Implementation of weighted Runge Kutta of order three

Hi,

I am trying to write my own code to implement the following Runge Kutta method:

The coefficients are given by:

I am using it to solve van Der Pol's Differential Equation (which is only dependant on y, hence I have not included the dependance on time for the coefficents:

I start by rewriting it as a first order system:

which gives:

My implementation looks like this:

clear all 
close all
clc

n = [320, 160, 80, 40, 20, 10]; % Amount of steps

egrid = zeros(1, 6); %For storage of the error for a given amount of steps
hgrid = zeros(1, 6); %For storage of the step size
y_NMax = 5.0087e-01; % Numerical solution with N = 320
counter = 0;
for N = n
    counter = counter + 1;
    T = 1;
    h = T/N; 
    t = linspace(0,T, N);
    u = zeros(2,length(t));
    u(1,1) = 1;
    u(2,1) = 0;

    for i = 1:(length(t) - 1)
        kvector = k(u(:,i), h);
        u(:,i + 1) = u(:,i) + h/6*(kvector);
    end
    
    hgrid(1, counter) = h;
    egrid(1, counter) = abs(u(1, end) - y_NMax);
    
    %plot(t, u(1,:));
    %hold on
    %plot(t,u(2,:));

end

% Code to plot and estimate the error

%loglog(hgrid(2:6), egrid(2:6))
%hold on
%loglog(hgrid(2:6), egrid(2:6), 'x')
%hold on
%loglog(hgrid(2:6), hgrid(2:6))
%hold on
%loglog(hgrid(2:6), hgrid(2:6).^2)
%hold on
%loglog(hgrid(2:6), hgrid(2:6).^3)


function dudt = uprime(u)
    du1dt = u(2);
    du2dt = -u(1) + u(2)*(1-u(1)^2);
    dudt = [du1dt; du2dt];
end

function kvector = k(u, h)
    k1 = uprime(u);
    k2 = uprime(u + h*k1);
    k3 = uprime(u + h/4*(k1 + k2));
    kvector = k1 + k2 + 4*k3 ;
end

I have tried to use the built in ODE45 to check my work, and it seems like I get a decent solution.

However, when I plot the error as a function of the step size h in the loglog-plot and compare it to the plot of h = h, I am led to belive that I have some error in my code, since it seems the error has a straight linear relationship to the step size h. This would not be in line with the theoretical results.

[Asad (Mehrzad) Khoddam]

It may doesn't change the error but for using linspace, change your code to:

t = linspace(0,T, N+1);

N+1 is the total number of points including start and end points, But N is the number oof divisions

[Felix Darke]

Hi,

Thanks for your reply. That did make my result get even closer to the result from ODE45, which is good.

However, I am still struggling with the error.

I get to take on the following values:

The loglog-plot for the error looks like this (the red line is h=h for comparison). With this result, it looks like the error is almost zero, and then becoming linear.

The error is approximated by

For comparision, ODE45 gives me