Improved Euler's method compared to 4th order Runge - Kutta's method

Is the improved Euler's method as good as the 4th order Runge - Kutta's method ? Or is it the step h that makes it almost not distinguishable ?

%--- Exercise 1 - A
%--- Solving the diff. eqn. N'(t) = N(t) - c*N^2(t) with c = 0.5 using
%--- Euler's method, Runge - Kutta method, and predictor-corrector method, and comparing in graph.
%--- I have four initial values for N0.

clc, clear all, close all

tmax = 50; % input('Enter tmax = ');
t = [0:tmax];
c = 0.5;
h = 0.1;
f = @(t,N) N - c*N.^2;
N0 = [0.1 0.5 1.0 2.0];
L = length(t)-1;
N = zeros(1,L);

for i = 1:length(N0)

figure(1)
subplot(2,2,i), hold on
title({'N'' = N - 0.5N^2',sprintf('N_0 = %.1f', N0(i))}),xlabel('t'), ylabel('N(t)'), hold on

%---------- Euler's method ----------

for j = 1:L
        
        EulerN(1) = N0(i);
        EulerN(1,j+1) = EulerN(j) + h*f(':',EulerN(j));
        
end

        plot(t,EulerN), hold on
        
%---------- Improved Euler's method ----------
        
for j = 1:L
        
        ImpN(1) = N0(i);
        ImpN(1,j+1) = ImpN(j) + 0.5*h*( f ( t(j), ImpN(j) ) + f ( t(j+1), ImpN(j) + h * f ( t(j), ImpN(j) ) ) );
         
end

        plot(t,ImpN,'r.'), hold on

%---------- Runge - Kutta method ----------
        
for j = 1:tmax
    
        RKN(1) = N0(i);
        
        K1 = f(t(j),RKN(j));
        K2 = f(t(j) + h*0.5, RKN(j) + h*K1*0.5);
        K3 = f(t(j) + h*0.5, RKN(j) + h*K2*0.5);
        K4 = f(t(j) + h, RKN(j) + h*K3);
    
        RKN(j+1) = RKN(j) + h*((K1 + K4)/6 + (K2 + K3)/3);
        
end

        plot(t,RKN,'k-.')
        legend({'Simple Euler','Improved Euler','Runge - Kutta'},'Location','Southeast')
        
end

[Torsten]

If the equation is non-stiff and you use fixed time stepping, usually the order of the method indicates which one is "best".

If the solutions are almost identical, you are correct: probably the stepsize is already small enough that the solutions lie one upon the other.

Answer 1

[John D'Errico]

Is it as good? COMPARE IT TO THE ANALYTICAL SOLUTION! We have gone through exactly this before.

A big problem however, is you are using the WRONG step size. That is, if you have

t = [0:tmax];

then your step size is 1. ONE. Do you see that? If you do that, then h must be 1. NOT 0.1. I'll use a step size of 0.1 though, which means I need to set t properly.

Next, you set the loop for Runge-Kutta to go from 1 to tmax. But the other loops went from 1 to L.

c = 0.5;
h = 0.1;
syms n(t);
n_t = dsolve(diff(n,1) == n - c*n^2,n(0) == 0.1)

n_t = 

nfun = matlabFunction(n_t)

nfun = function_handle with value:
    @(t)(exp(t).*2.0)./(exp(t)+1.9e+1)

That is the analytical solution. I'll strip out some of what you have written and solve lt for one initial value.

tmax = 10;
t = [0:h:tmax];
f = @(t,N) N - c*N.^2;
N0 = 0.1;

L = length(t)-1;
for j = 1:L
        EulerN(1) = N0;
        EulerN(1,j+1) = EulerN(j) + h*f(':',EulerN(j));
end
        
%---------- Improved Euler's method ----------
for j = 1:L
        ImpN(1) = N0;
        ImpN(1,j+1) = ImpN(j) + 0.5*h*( f ( t(j), ImpN(j) ) + f ( t(j+1), ImpN(j) + h * f ( t(j), ImpN(j) ) ) );
end

%---------- Runge - Kutta method ----------        
for j = 1:L
        RKN(1) = N0;
        K1 = f(t(j),RKN(j));
        K2 = f(t(j) + h*0.5, RKN(j) + h*K1*0.5);
        K3 = f(t(j) + h*0.5, RKN(j) + h*K2*0.5);
        K4 = f(t(j) + h, RKN(j) + h*K3);
        
        RKN(j+1) = RKN(j) + h*((K1 + K4)/6 + (K2 + K3)/3);
end

N = nfun(t);

plot(t,N,'k-',t,EulerN,'g-',t,ImpN,'b-',t,RKN,'r-')
title 'Predictions'
grid on
legend({'True solution','Simple Euler','Improved Euler','Runge - Kutta'},'Location','Southeast')

plot(t,N - EulerN,'g-',t,N - ImpN,'b-',t,N - RKN,'r-')
title('Errors')
grid on
legend({'Simple Euler','Improved Euler','Runge - Kutta'},'Location','Northeast')

You should see that now, in the first plot, ALL of the curves pretty much lie on top of each other, now that I used comparable step sizes, etc. The problem is, you can't really see the difference. BUT in the error plot, how do they compare?

You should then see that simple Euler is ok, but it is way worse than improved Euler, and Runge-Kutta pretty much nails it, and does so far better than did either version of Euler. You really need to compare the methods using the same step sizes though. Else you will get meaningless, uninterpretable garbage.

You really need to learn to be more careful in your code. Otherwise, how can you trust the results you get?

[Left Terry]

@John D'Errico What does the line "nfun = matlabFunction(n_t)" do ? Also what "N = nfun(t)" does ?

[John D'Errico]

matlabFunction takes a symbolic expression (in this case, the exact solution in n_t), and converts it to a function you can evaluate.

So now the line:

N = nfun(t);

evaluates the function I just created, at the list of points in the vector t.

[Left Terry]

@John D'Errico I got it, thank you. It almost works in my code, except from the last graph for the errors.

%--- Exercise 1 - A
%--- Solving the diff. eqn. N'(t) = N(t) - c*N^2(t) with c = 0.5 using
%--- Euler's method, Runge - Kutta method, and predictor-corrector method, and comparing in graph.
%--- I have four initial values for N0.

clc, clear all, close all

tmax = 10; % input('Enter tmax = ');
c = 0.5;
h = 0.1;
f = @(t,N) N - c*N.^2;
N0 = [0.1 0.5 1.0 2.0];

for i = 1:length(N0)
    
syms n(t)
Df = diff(n) == n - c*n^2;
nAnalytical = dsolve(Df,n(0) == N0(i));

nfun = matlabFunction(nAnalytical);

figure(i)
subplot(1,2,1)
fplot(nAnalytical,[0 tmax],'k.'), hold on
title({'N'' = N - 0.5N^2',sprintf('N_0 = %.1f', N0(i))}),xlabel('t'), ylabel('N(t)'), hold on

t = [0:h:tmax];
L = length(t)-1;
N = zeros(1,L);

%---------- Simple Euler's method ----------

for j = 1:L
        
        EulerN(1) = N0(i);
        EulerN(1,j+1) = EulerN(j) + h*f(':',EulerN(j));
               
end

        plot(t,EulerN), hold on
        
%---------- Improved Euler's method ----------
        
for j = 1:L
        
        ImpN(1) = N0(i);
        ImpN(1,j+1) = ImpN(j) + 0.5*h*( f ( t(j), ImpN(j) ) + f ( t(j+1), ImpN(j) + h * f ( t(j), ImpN(j) ) ) );
         
end

        plot(t,ImpN), hold on

%---------- Runge - Kutta method ----------
        
for j = 1:L
    
        RKN(1) = N0(i);
        
        K1 = f(t(j),RKN(j));
        K2 = f(t(j) + h*0.5, RKN(j) + h*K1*0.5);
        K3 = f(t(j) + h*0.5, RKN(j) + h*K2*0.5);
        K4 = f(t(j) + h, RKN(j) + h*K3);
    
        RKN(j+1) = RKN(j) + h*((K1 + K4)/6 + (K2 + K3)/3);
        
end
       plot(t,RKN)
       legend({'Analytical solution','Simple Euler','Improved Euler','RK4'},'Location','Southeast')
       
       N = nfun(t);
       subplot(1,2,2), plot(t,abs(N - EulerN)*100,'k-',t,abs(N - ImpN)*100,'b-',t,abs(N - RKN)*100,'r-')
       title('Error of each method'), ylabel('Error (%)'), legend({'Simple Euler','Improved Euler','RK4'})
end

Error using sym/matlabFunction>@()2.0
Too many input arguments.

[Walter Roberson]

tmax = 10; % input('Enter tmax = ');
c = 0.5;
h = 0.1;
f = @(t,N) N - c*N.^2;
N0 = [0.1 0.5 1.0 2.0];
for i = 1:length(N0)
    N0(i)
    syms n(t)
    Df = diff(n) == n - c*n^2;
    nAnalytical = dsolve(Df,n(0) == N0(i))
end

ans = 0.1000

nAnalytical = 

ans = 0.5000

nAnalytical = 

ans = 1

nAnalytical = 

ans = 2

nAnalytical = 

2

we can see that when N0 becomes 2, the analytic solution is the constant function n(t) = 2. But that constant function does not involve any symbolic variables, so when you

nfun = matlabFunction(nAnalytical)

you get out @() 2 as the function. That anonymous function does not accept any input parameters.

To get this to work, you need to use

nfun = matlabFunction(nAnalytical, 'vars', t)

[Left Terry]

@Walter Roberson Works perfectly. Thank you.