Effect of Step Size ODE

Hello,

my big deal is to demonstrate the effect of the step size of differential equations with an intial value problem.

• Make a choice for two methods of different orders (like Euler combined with Hern), (in total we have 4 Methods - Euler, Heun, Runge-Kutta, Simpson)
• start to apply both methods in order to compute (t_k,y_k) from (t_(k-1),y_(k-1)) with both methods using a step size h,
• compute the difference between the two values y_k (Euler) and y_k (Heun); this is called local error e_k here,
• When the local error is high (let us say, e_k> 0.01), the step size has to be reduced (for example h is replecd by h/2), and the computation is redone. Otherwise, when the local error is very small (let us say, e_k > 0.0001), the step size is increased (for example h is replaced by 1.5*h), and the computation of the y_k is redone with the new step size.

1. how can i combine the user input, if the user choose two variables? --> I think the switch case is not correct....
2. i have an error in the functions, matlab told me, that y has too many output arguments and i should consider the preallocating, but in the beginning, i defined y=zeros; what can i do?

clear all;

clc;

clf;

F=@(t,y) 200*(y+0.01).^2*(0.2-y); %Initial Value Problem

m = input ('Please choose the first method (1=Euler; 2=classical Runge-Kutta; 3=modified Simpson; 4=Heun: ');

l = input ('Please choose the second method (1=Euler; 2=classical Runge-Kutta; 3=modified Simpson; 4=Heun: ');

a = input('Enter left end ponit, a: ');

b = input('Enter right end point, b: ');

N = input('Enter no. of subintervals, n: ');

y0 = input('Enter the initial value problem, y0: ');

h = (b-a)/N; %Stepsize

x=a:h:b;

y=zeros;

% User Choice

switch m

case 1

plot(x,euler,'--','red')

case 2

plot(x,runge,'--','green')

case 3

plot(x,simpson,'--','black')

case 4

plot(x,heun,'--','blue')

end

switch l

case 1

plot(x,euler,'--','red')

case 2

plot(x,runge,'--','green')

case 3

plot(x,simpson,'--','black')

case 4

plot(x,heun,'--','blue')

end

% Euler´s Method

function euler

for i=1:(lenght(x)-1)

y(i+1)=y(i)+h*F(x(i),y(i));

end

end

% Runge-Kutta Method

function runge

for i=1:(lenght(x)-1)

k1 = F(x(i-1),y(i-1));

k2 = F(x(i-1)+0.5*h,y(i-1)+0.5*h*k1);

k3 = F((x(i-1)+0.5*h),(y(i-1)+0.5*h*k2));

k4 = F((x(i-1)+h),(y(i-1)+k3*h));

y(i) = y(i-1) + (h/6)*(k1+2*k2+2*k3+k4);

end

end

% modified Simpson Method

function simpson

for i=1:(lenght(x)-1)

k1=F(x(i-1),y(i-1));

k2=F(x(i-1)+(h/2),y(i-1)+(h/2)*k1);

k3=F(x(i-1)+h,y(i-1)+h*k2);

y(i)=y(i-1)+(h/6)*(k1+4*k2+k3);

end

end

% Heuns's Method

function heun

for i=1:(lenght(x)-1)

k1=F(x(i-1),y(i-1));

k2=F(x(i-1)+h,y(i-1)+h*k1);

y(i)=y(i-1)+(h/2)*(k1+k2);

end

end

Answer 1

[Alan Stevens]

Quite a few changes needed!

F=@(t,y) 200*(y+0.01).^2*(0.2-y); %Initial Value Problem
 
m = input ('Please choose the first method (1=Euler; 2=classical Runge-Kutta; 3=modified Simpson; 4=Heun: ');
l = input ('Please choose the second method (1=Euler; 2=classical Runge-Kutta; 3=modified Simpson; 4=Heun: ');
a = input('Enter left end point, a:  ');
b = input('Enter right end point, b:  ');
N = input('Enter no. of subintervals, n: ');                
y0 = input('Enter the initial value, y0:  ');
 
h = (b-a)/N; %Stepsize                   
x=a:h:b;
y=zeros(1,length(x));  y(1) = y0;

figure
hold on
% User Choice
switch m
    case 1
        ye = euler(x, y, F, h);
        plot(x,ye,'r--')
    case 2
        yr = runge(x, y, F, h);
        plot(x,yr,'g--')
    case 3
        ys = simpson(x, y, F, h);
        plot(x,ys,'k--')
    case 4
        yh = heun(x, y, F, h);
        plot(x,yh,'b--')
end  
 
switch l
    case 1
        ye = euler(x, y, F, h);
        plot(x,ye,'r--')
    case 2
        yr = runge(x, y, F, h);
        plot(x,yr,'g--')
    case 3
        ys = simpson(x, y, F, h);
        plot(x,ys,'k--')
    case 4
        yh = heun(x, y, F, h);
        plot(x,yh,'b--')
end  
hold off    
 
 
% Euler´s Method
function y = euler(x, y, F, h) 
        for i=2:(length(x))
            y(i)=y(i-1)+h*F(x(i-1),y(i-1));
        end
end
    
% Runge-Kutta Method
     
function y = runge(x, y, F, h)
        for i=2:(length(x))
         k1 = F(x(i-1),y(i-1));
         k2 = F(x(i-1)+0.5*h,y(i-1)+0.5*h*k1);
         k3 = F((x(i-1)+0.5*h),(y(i-1)+0.5*h*k2));
         k4 = F((x(i-1)+h),(y(i-1)+k3*h));
         y(i) = y(i-1) + (h/6)*(k1+2*k2+2*k3+k4);
        end
end
     
 
% modified Simpson Method
function y = simpson(x, y, F, h)
     for i=2:(length(x))
         k1=F(x(i-1),y(i-1));
         k2=F(x(i-1)+(h/2),y(i-1)+(h/2)*k1);
         k3=F(x(i-1)+h,y(i-1)+h*k2);
         y(i)=y(i-1)+(h/6)*(k1+4*k2+k3);
     end
end     
 
% Heuns's Method 
function y = heun(x, y, F, h)
     for i=2:(length(x))
         k1=F(x(i-1),y(i-1));
         k2=F(x(i-1)+h, y(i-1)+h*k1);
         y(i)=y(i-1)+(h/2)*(k1+k2);
     end
end