Runge kutta method Matlab 5th order

Question

How can we moditify the code of Runge Kutta method 2th order from site http://www.computeraideddesignguide.com/euler-vs-runge-kutta-circle-matlab-code/ into the code of Runge Kutta method 5th order?

i tried, but returns back like a spiral. i used Dormand Prince tableau. Can you help to find my mistakes in my code?

 clc;
 clear all;
 h=0.01;
 N=5000;
 x=zeros(1,N);
 y=zeros(1,N);
 x(1)=1; %initial value x(0)
 y(1)=0; %initial value y(0)
 %iteration with N step
 % Forward Euler method
 c2=0.2; c3 = 0.3 ;c4=0.8 ;c5=8./9.; c6=1.; c7=1.;
 a21 = 0.2; a31 = 3./40.0; a32 = 9.0/40.0; a41 = 44.0/45.0;  a42=56.0/15.0; a43=32.0/9.0;
 a51=19372.0/6561.0; a52=25360.0/2187.0; a53=64448.0/6561.0; a54=212.0/729.0;
 a61=9017.0/3168.; a62=355.0/33.0; a63=46732.0/5247.0; a64=49.0/176.0; a65=5103.0/18656.0;
 a71=35.0/384.0; a73=500.0/1113.0; a74=125.0/192.0 ;a75=2187.0/6784.0; a76=11.0/84.0;
 for i=1:N
 xk1=y(i);
 yk1=-x(i);
 xk2 = y(i) + c2*h*(a21*xk1);
 yk2 = -(x(i) + c2*h*(a21*yk1));
 xk3 = y(i) + c3*h*(a31*xk1 + a32*xk2);
 yk3 = -(x(i) + c3*h*(a31*yk1 + a32*yk2));
 xk4 = y(i) + c4*h*(a41*xk1 - a42*xk2 + a43*xk3);
 yk4 = -(x(i) + c4*h*(a41*yk1 - a42*yk2 + a43*yk3));
 xk5 = y(i) + c5*h*(a51*xk1 - a52*xk2 + a53*xk3 - a54*xk4);
 yk5 = -(x(i) + c5*h*(a51*yk1 - a52*yk2 + a53*yk3 - a54*yk4 ));
 xk6 = y(i) + c6*h*(a61*xk1 - a62*xk2 + a63*xk3 + a64*xk4 - a65*xk5);
 yk6 = -(x(i) + c6*h*(a61*yk1 - a62*yk2 + a63*yk3 + a64*yk4 - a65*yk5));
 xk7 = y(i) + c7*h*(a71*xk1 + a73*xk3 + a74*xk4 - a75*xk5 + a76*xk6);
 yk7 = -(x(i) + c7*h*(a71*yk1 + a73*yk3 + a74*yk4 - a75*yk5 + a76*yk6));
 x(i+1) = x(i) + h *((5179.0/57600.)*xk1 + (7571.0/16695.0)*xk3 + (393.0/640.0)*xk4 - (92097.0/339200.0)*xk5 +(187.0/2100.0)*xk6 + (1.0/40)*xk7);
 y(i+1) = y(i) + h * ((5179.0/57600.)*yk1 + (7571.0/16695.0)*yk3 + (393.0/640.0)*yk4 - (92097.0/339200.0)*yk5 +(187.0/2100.0)*yk6 + (1.0/40)*yk7);
 end
 angle=linspace(0,2*pi); % used to plot the unit circle
 plot(x,y,'r',sin(angle),cos(angle),'b')
 legend('Range Kutta', 'Unit circle')

Answer 1

You need to use the xk1,xk2,... as updates for the x and the yk1,yk2,... as updates for the y. At the moment you have that mixed up, which reduces the order of the method to 1 and you get spirals like with the Euler method. Also, you are using the c values wrong. As the ODE does not depend on time, they should not even be present.

 xk2 = y(i) + h*(a21*yk1);
 yk2 = -(x(i) + h*(a21*xk1));
 xk3 = y(i) + h*(a31*yk1 + a32*yk2);
 yk3 = -(x(i) + h*(a31*xk1 + a32*xk2));

etc.

You will find that, with a correct implementation, you can dramatically increase the step size to only a few steps per circle, about 10 or less. In that case you would need the "dense output" interpolation for additional values per step.