Implementing Runge-Kutta 4

Question

I am trying to implement Runge-Kutta 4 to solve a differential equation and need some insight.

The main problem that I'm having is that my error for these values is around 0.09, when it should be around 0.0001*10-4 , so there's something wrong with my implementation of rk4, but I have no idea where I'm making the error. If you could point me in the direction of where my implementation of runge-kutta is off, that would be great. (Note, we are able to compute the error because the solution is a circle, so I know that the end point should be the same as the initial condition, (1,0), and the error is the distance between the computed endpoint and (1,0). This assignment is for learning numerical methods.

I repeat: I AM NOT LOOKING FOR A SOLUTION, I am looking for someone to help point me in the direction of what's going wrong in my code, because I've been staring at this for god knows how long...

How the code works: between the function declaration and for loop is setting the initial values (again, p and q are irrelevant in this part of the problem), the for loop is iterating and solving the function numerically. I use runge kutta 4 as described on the wikipedia article.

The code I wrote is written is below and N=400 (400 steps), T=42 (terminal time of 4pi2), (x,y)=(1,0) (initial conditions), and gamma = 1 (gamma is a parameter in the equation), and (p,q) is irrelevant for the part about which I am asking. P1PC is the .m file that contains the differential equation, and is also below.

function rk(N,T,x,y,gamma,p,q)
timestep=T/N;
xy=zeros(N,4);
xy(1,:)=[x,y,p,q];
k=zeros(4,2);%k(:,1) is for x, k(:,2) is for y
for index=2:N
    [k(1,1),k(1,2)]=P1PC(gamma,xy(index-1,1),xy(index-1,2));
    [k(2,1),k(2,2)]=P1PC(gamma,xy(index-1,1)+0.5*k(1,1)*timestep,xy(index-1,2)+0.5*k(1,2)*timestep);
    [k(3,1),k(3,2)]=P1PC(gamma,xy(index-1,1)+0.5*k(2,1)*timestep,xy(index-1,2)+0.5*k(2,2)*timestep);
    [k(4,1),k(4,2)]=P1PC(gamma,xy(index-1,1)+k(3,1)*timestep,xy(index-1,2)+k(3,2)*timestep);
    xy(index,1)=xy(index-1,1)+(timestep/6)*(k(1,1)+2*k(2,1)+2*k(3,1)+k(4,1));
    xy(index,2)=xy(index-1,2)+(timestep/6)*(k(1,2)+2*k(2,2)+2*k(3,2)+k(4,2));
end
plot(xy(:,1),xy(:,2));
error=sqrt((1-xy(N,1))^2+(0-xy(N,2))^2)
xy(N,1)
xy(N,2)
end

Here is the matlab code for function I am solving (P1PC):

function [kx,ky]=P1PC(gamma,x,y)
    r=x^2+y^2;
    ky=(gamma*x)/(2*pi*r);
    kx=(gamma*(-y))/(2*pi*r);
    end

yes, the kx and ky are my shorthand for dx/dt, dy/dt, although i supposed it would make more sense to write dx and dy lol — malxmusician212
Check your end time. You are taking N-1 steps of h=T/N each. I think you need for index=2:N+1 to get to the final values. — John Alexiou
Your end time does not correspond to getting back to [x,y]=[1,0]. To complete 2*pi radians of rotation, T must be (2*pi)^2 which is roughly 39.478, not 42. — David

John Alexiou John Alexiou · Accepted Answer · 2015-01-20T23:25:49

I see two issues. One is that the end time is only reached after taking N steps, and your code takes N-1 steps. Most importantly, your definition of error is wrong. You want to check if the radius is one and hence error=sqrt(x^2+y^2)-1

See the code below I used (a bit simplified) to test the algorithm

P1PC =@(gamma,x,y)[-gamma*y,gamma*x]/(2*pi*(x^2+y^2));
T = 42;
N = 400;
h=T/N;
time=0;
x0=1;
y0=0;
gamma=1;
xy = zeros(N+1,2);
xy(1,:) = [x0,y0];
for i=2:N+1
    k1 = P1PC(gamma, xy(i-1,1),xy(i-1,2));
    k2 = P1PC(gamma, xy(i-1,1)+(h/2)*k1(1), xy(i-1,2)+(h/2)*k1(2));
    k3 = P1PC(gamma, xy(i-1,1)+(h/2)*k2(1), xy(i-1,2)+(h/2)*k2(2));
    k4 = P1PC(gamma, xy(i-1,1)+(h)*k3(1), xy(i-1,2)+(h)*k3(2));
    xy(i,:) = xy(i-1,:) + (h/6)*(k1+2*k2+2*k3+k4);
    time = time + h;
end

plot(xy(:,1),xy(:,2));
disp(['time=',num2str(time)])
error=sqrt(xy(N+1,1)^2+xy(N+1,2)^2)-1;
disp(['error=',num2str(error)])

Produces the output

time=42
error=3.0267e-011

Implementing Runge-Kutta 4

1 Answers