4th order Runge - Kutta vs. 4th order Predictor - Corrector

Hello once again. As we can see, the error of the 4th order predictor - corrector (assuminig my code for this is correct) is greater than that of 4th order Runge - Kutta. Shouldn't be the other way around ?

%--- 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, close all, format long

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];
LN = length(N0);
%---------- Analytical Solution ----------

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

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

%-------------------------------------------

t = 0:h:tmax;
L = length(t)-1;

%---------- Simple Euler's method ----------
    
    EulerN = zeros(1,L);
    
    for j = 1:L
        
        EulerN(1) = N0(i);
        EulerN(1,j+1) = EulerN(j) + h*f(':',EulerN(j));
                      
    end

% ---------- Improved Euler's method ----------
        
    ImpN = zeros(1,L);

    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

%---------- Runge - Kutta method ----------
        
    RKN = zeros(1,L);

    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/6*(K1 + K4 + 2*K2 + 2*K3);
        
    end
    
%---------- 4th order Predictor - Corrector ----------

    PredN = zeros(1,L);
    CorN = zeros(1,L);
    for k = 1:L+1
       if k <= 4
        PredN(k) =  RKN(k);
        CorN(k) =  RKN(k);
        else
        PredN(k) = PredN(k-1) + h/24*( 55*f(t(k-1),RKN(k-1)) - 59*f(t(k-2),RKN(k-2)) + 37*f(t(k-3),RKN(k-3)) - 9*f(t(k-4),RKN(k-4)));
        CorN(k) = PredN(k-1) + h/24*(9*f(t(k),PredN(k)) + 19*f(t(k-1),RKN(k-1)) - 5*f(t(k-2),RKN(k-2)) + f(t(k-3),RKN(k-3)));
       end
    end
    PC4 = CorN;
%---------------------------------------------------
%---------- Plotting -------------------------------
       
       figure(i)
       
       subplot(3,1,1)
       fplot(nAnalytical,[0 tmax],'k')
       title({'Solution of N'' = N - 0.5N^2',sprintf('for the initial condition N_0 = %.1f', N0(i))})
       xlabel('t'), ylabel('N(t)'), hold on, grid on

       plot(t,EulerN,t,ImpN,t,RKN,t,PC4)
       legend({'Analytical solution','Simple Euler','Improved Euler','RK4','PC4'},'Location','Southeast'), grid on
       
       N = nfun(t);
       
       subplot(3,1,2)
       plot(t,abs(N - EulerN),t,abs(N - ImpN),t,abs(N - RKN),t,abs(N - PC4))
       title(sprintf('Error of each method for N_0 = %.1f',N0(i)))
       xlabel('t'), ylabel('Error')
       legend({'Simple Euler','Improved Euler','RK4','PC4'}), grid on, hold on
       
       subplot(3,1,3)
       plot(t,abs(N - RKN),t,abs(N - PC4))
       title(sprintf('Error of RK4 and PC4 for N_0 = %.1f',N0(i)))
       xlabel('t'), ylabel('Error')
       legend({'RK4','PC4'}), grid on, hold on

%----------------------------------------------------------
end

Answer 1

[Torsten]

%--- 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, format long
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];
LN = length(N0);
%---------- Analytical Solution ----------
for i = 1:LN
    
syms n(t)
Df = diff(n) == n - c*n^2;
nAnalytical = dsolve(Df,n(0) == N0(i));
nfun = matlabFunction(nAnalytical, 'vars', t);
% %-------------------------------------------
t = 0:h:tmax;
L = length(t)-1;
%---------- Simple Euler's method ----------
    
    EulerN = zeros(1,L);
    
    for j = 1:L
        
        EulerN(1) = N0(i);
        EulerN(1,j+1) = EulerN(j) + h*f(':',EulerN(j));
                      
    end
% ---------- Improved Euler's method ----------
        
    ImpN = zeros(1,L);
    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
%---------- Runge - Kutta method ----------
        
    RKN = zeros(1,L);
    RKN(1) = N0(i);
    
    for j = 1:L
           
        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/6*(K1 + K4 + 2*K2 + 2*K3);
        
    end
    
%---------- 4th order Predictor - Corrector ----------
    PC4 = zeros(1,L+1);
    PC4(1:4) = RKN(1:4);
    for j = 4:L
       PC40 = PC4(j)+h*(55.0*f(t(j),PC4(j))-59.0*f(t(j-1),PC4(j-1))+37.0*f(t(j-2),PC4(j-2))-9.0*f(t(j-3),PC4(j-3)))/24.0; 
       % correct PC4(j+1)
       PC4(j+1) = PC4(j)+h*(9.0*f(t(j+1),PC40)+19.0*f(t(j),PC4(j))-5.0*f(t(j-1),PC4(j-1))+f(t(j-2),PC4(j-2)))/24.0; 
    end
     
%---------- Plotting -------------------------------
       
       figure(i)
       
       subplot(3,1,1)
       fplot(nAnalytical,[0 tmax],'k')
       title({'Solution of N'' = N - 0.5N^2',sprintf('N_0 = %.1f', N0(i))})
       xlabel('t'), ylabel('N(t)'), hold on, grid on
       plot(t,EulerN,t,ImpN,t,RKN,t,PC4)
       legend({'Analytical solution','Simple Euler','Improved Euler','RK4','PC4'},'Location','Southeast'), grid on
       
       N = nfun(t);
       
       subplot(3,1,2)
       plot(t,abs(N - EulerN),t,abs(N - ImpN),t,abs(N - RKN),t,abs(N - PC4))
       title(sprintf('Error of each method for N_0 = %.1f',N0(i)))
       xlabel('t'), ylabel('Error')
       legend({'Simple Euler','Improved Euler','RK4','PC4'}), grid on, hold on
       
       subplot(3,1,3)
       plot(t,abs(N - RKN),t,abs(N - PC4))
       title(sprintf('Error of RK4 and PC4 for N_0 = %.1f',N0(i)))
       xlabel('t'), ylabel('Error')
       legend({'RK4','PC4'}), grid on, hold on
%----------------------------------------------------------
end

●

[Left Terry]

@Torsten So, PC4 it's still worst than RK4, but is this true in general ?

[Torsten]

So, PC4 it's still worst than RK4, but is this true in general ?

Not in all regions for t. And both have order 4 - so why do you think PC4 should be more precise than RK4 ?

[Torsten]

In order to avoid unnecessary evaluations of the derivative function, I suggest using the following code for the predictor-corrector scheme:

%---------- 4th order Predictor - Corrector ----------
    PC4 = zeros(1,L+1);
    PC4(1:4) = RKN(1:4);
    fm1 = f(t(3),PC4(3));
    fm2 = f(t(2),PC4(2));
    fm3 = f(t(1),PC4(1));
    for j = 4:L
       fm0 = f(t(j),PC4(j));  
       PC40 = PC4(j)+h*(55.0*fm0-59.0*fm1+37.0*fm2-9.0*fm3)/24.0; 
       % correct PC4(j+1)
       fPC40 = f(t(j+1),PC40);
       PC4(j+1) = PC4(j)+h*(9.0*fPC40+19.0*fm0-5.0*fm1+fm2)/24.0; 
       fm3 = fm2;
       fm2 = fm1;
       fm1 = fm0;
    end

Answer 2

[Torsten]

For k>4, the variable RKN should no longer appear in the equations

        PredN(k) = PredN(k-1) + h/24*( 55*f(t(k-1),RKN(k-1)) - 59*f(t(k-2),RKN(k-2)) + 37*f(t(k-3),RKN(k-3)) - 9*f(t(k-4),RKN(k-4)));
        CorN(k) = PredN(k-1) + h/24*(9*f(t(k),PredN(k)) + 19*f(t(k-1),RKN(k-1)) - 5*f(t(k-2),RKN(k-2)) + f(t(k-3),RKN(k-3)));

[Left Terry]

@Torsten Correct. I think then that the correct formulae is the following

PredN(k) = CorN(k-1) + h/24*( 55*f(t(k-1),CorN(k-1)) - 59*f(t(k-2),CorN(k-2)) + 37*f(t(k-3),CorN(k-3)) - 9*f(t(k-4),CorN(k-4)));
CorN(k) = PredN(k-1) + h/24*(9*f(t(k),PredN(k)) + 19*f(t(k-1),PredN(k-1)) - 5*f(t(k-2),PredN(k-2)) + f(t(k-3),PredN(k-3)));

[Left Terry]

@Torsten No, the above formulas doesn't fix the problem.

%--- 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, format long

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];
LN = length(N0);

%---------- Analytical Solution ----------

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

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

% %-------------------------------------------

t = 0:h:tmax;
L = length(t)-1;

%---------- Simple Euler's method ----------
    
    EulerN = zeros(1,L);
    
    for j = 1:L
        
        EulerN(1) = N0(i);
        EulerN(1,j+1) = EulerN(j) + h*f(':',EulerN(j));
                      
    end

% ---------- Improved Euler's method ----------
        
    ImpN = zeros(1,L);

    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

%---------- Runge - Kutta method ----------
        
    RKN = zeros(1,L);

    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/6*(K1 + K4 + 2*K2 + 2*K3);
        
    end
    
%---------- 4th order Predictor - Corrector ----------

    PredN = zeros(1,L);
    CorN = zeros(1,L);
    for k = 1:L+1
       if k <= 4
        PredN(k) =  RKN(k);
        CorN(k) =  RKN(k);
        else
        PredN(k) = CorN(k-1) + h/24*( 55*f(t(k-1),CorN(k-1)) - 59*f(t(k-2),CorN(k-2)) + 37*f(t(k-3),CorN(k-3)) - 9*f(t(k-4),CorN(k-4)));
        CorN(k) = PredN(k-1) + h/24*(9*f(t(k),PredN(k)) + 19*f(t(k-1),PredN(k-1)) - 5*f(t(k-2),PredN(k-2)) + f(t(k-3),PredN(k-3)));
       end
       
%         fprintf('\n%3d)\tRunge - Kutta\t\tPrediction\t\tCorrection\n',k);
%         fprintf('\t\t%.8f\t\t\t%.8f\t\t%.8f\n',RKN(k),PredN(k),CorN(k));
    end
    PC4 = CorN;
    
%---------------------------------------------------
%---------- Plotting -------------------------------
       
       figure(i)
       
       subplot(3,1,1)
       fplot(nAnalytical,[0 tmax],'k')
       title({'Solution of N'' = N - 0.5N^2',sprintf('N_0 = %.1f', N0(i))})
       xlabel('t'), ylabel('N(t)'), hold on, grid on

       plot(t,EulerN,t,ImpN,t,RKN,t,PC4)
       legend({'Analytical solution','Simple Euler','Improved Euler','RK4','PC4'},'Location','Southeast'), grid on
       
       N = nfun(t);
       
       subplot(3,1,2)
       plot(t,abs(N - EulerN),t,abs(N - ImpN),t,abs(N - RKN),t,abs(N - PC4))
       title(sprintf('Error of each method for N_0 = %.1f',N0(i)))
       xlabel('t'), ylabel('Error')
       legend({'Simple Euler','Improved Euler','RK4','PC4'}), grid on, hold on
       
       subplot(3,1,3)
       plot(t,abs(N - RKN),t,abs(N - PC4))
       title(sprintf('Error of RK4 and PC4 for N_0 = %.1f',N0(i)))
       xlabel('t'), ylabel('Error')
       legend({'RK4','PC4'}), grid on, hold on

%----------------------------------------------------------
end

%--- Exercise 1 - B
%--- Solving the diff. eqn. y'(t) = y(t) - c*y^2(t) with initial condition
%--- various values for c using Runge - Kutta method

% clc
% c = [0.1 0.2 0.5 1.0 2.0 5.0];
% Nmax = 1./c;
% 
% % t(1) = 0;
% L = length(t);
% N = zeros(1,L);
% N(1) = 0.1;
% 
% for i = 1:length(c)
% 
%     f = @(t,N) N - c(i)*N^2;
%     
%     for j = 1:L-1
%         
%         K1 = f(t(j),N(j));
%         K2 = f(t(j) + h*0.5, N(j) + h*K1*0.5);
%         K3 = f(t(j) + h*0.5, N(j) + h*K2*0.5);
%         K4 = f(t(j) + h, N(j) + h*K3);
%     
%         N(j+1) = N(j) + h/6*(K1 + K4 + 2*K2 + 2*K3); 
%     end
%     figure(LN+1)
%     subplot(2,1,1)
%     plot(t,N), hold on, grid on
%     title('Solutions of N'' = N - cN^2 for different values of c'), xlabel('t'), ylabel('N(t)')
% end
% legend({'c = 0.1','c = 0.2','c = 0.5','c = 1','c = 2','c = 5'},'Location','Best')
% subplot(2,1,2)
% plot(c,Nmax), title('N_{max} vs. c'), xlabel('c'), ylabel('N_{max}'), grid on, text(0.5,4,'N_{max} = 1/c')