Ndtyjky

I´m attempting to use Runge-Kutta method to compare it to the lsode function. But it is performing rather poorly, every other method I used (fowards and backwards Euler, Heun) to compare to lsode do a way better job to the point they are almost indistinguishable from lsode.

This is what my code returns
https://i.stack.imgur.com/vJ6Yi.png

If anyone can pointout a way to improve it or if I doing something wrong I´d appreciate it.

The following is what I use for the Runge-Kutta method

%Initial conditions



u(1) = 1;

v(1) = 2;

p(1) = -1/sqrt(3);

q(1) = 1/sqrt(3);



%Graf interval / step size

s0 = 0;

sf = 50;

h = 0.25;



n=(sf-s0)/h;



s(1) = s0;



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



for j = 2:n



  i = j-1;



  k1_u(j) = p(i);

  k1_v(j) = q(i);

  k1_p(j) = (-2*v(i)*p(i)*q(i)) / (u(i)*u(i) + v(i)*v(i) + 1);

  k1_q(j) = (-2*u(i)*p(i)*q(i)) / (u(i)*u(i) + v(i)*v(i) + 1);



  u1(j) = p(i) + (1/2)*k1_u(j)*h;

  v1(j) = q(i) + (1/2)*k1_v(j)*h;

  p1(j) = (-2*v(i)*p(i)*q(i)) / (u(i)*u(i) + v(i)*v(i) + 1) + (1/2)*k1_p(j)*h;

  q1(j) = (-2*u(i)*p(i)*q(i)) / (u(i)*u(i) + v(i)*v(i) + 1) + (1/2)*k1_q(j)*h;



  k2_u(j) = p1(j);

  k2_v(j) = q1(j);

  k2_p(j) = (-2*v1(j)*p1(j)*q1(j)) / (u1(j)*u1(j) + v1(j)*v1(j) + 1);

  k2_q(j) = (-2*u1(j)*p1(j)*q1(j)) / (u1(j)*u1(j) + v1(j)*v1(j) + 1);



  u2(j) = p(i) + (1/2)*k2_u(j)*h;

  v2(j) = q(i) + (1/2)*k2_v(j)*h;

  p2(j) = (-2*v(i)*p(i)*q(i)) / (u(i)*u(i) + v(i)*v(i) + 1) + (1/2)*k2_p(j)*h;

  q2(j) = (-2*u(i)*p(i)*q(i)) / (u(i)*u(i) + v(i)*v(i) + 1) + (1/2)*k2_q(j)*h;



  k3_u(j) = p2(j);

  k3_v(j) = q2(j);

  k3_p(j) = (-2*v2(j)*p2(j)*q2(j)) / (u2(j)*u2(j) + v2(j)*v2(j) + 1);

  k3_q(j) = (-2*u2(j)*p2(j)*q2(j)) / (u2(j)*u2(j) + v2(j)*v2(j) + 1);



  u3(j) = p(i) + k3_u(j)*h;

  v3(j) = q(i) + k3_v(j)*h;

  p3(j) = (-2*v(i)*p(i)*q(i)) / (u(i)*u(i) + v(i)*v(i) + 1) + k3_p(j)*h;

  q3(j) = (-2*u(i)*p(i)*q(i)) / (u(i)*u(i) + v(i)*v(i) + 1) + k3_q(j)*h;



  k4_u(j) = p3(j);

  k4_v(j) = q3(j);

  k4_p(j) = (-2*v3(j)*p3(j)*q3(j)) / (u3(j)*u3(j) + v3(j)*v3(j) + 1);

  k4_q(j) = (-2*u3(j)*p3(j)*q3(j)) / (u3(j)*u3(j) + v3(j)*v3(j) + 1);





    s(j) = s(j-1) + h;

    u(j) = u(j-1) + (h/6)*(k1_u(j) + 2*k2_u(j) + 2*k3_u(j) + k4_u(j));

    v(j) = v(j-1) + (h/6)*(k1_v(j) + 2*k2_v(j) + 2*k3_v(j) + k4_v(j));

    p(j) = p(j-1) + (h/6)*(k1_p(j) + 2*k2_p(j) + 2*k3_p(j) + k4_p(j));

    q(j) = q(j-1) + (h/6)*(k1_q(j) + 2*k2_q(j) + 2*k3_q(j) + k4_q(j));



endfor



subplot(2,3,1), plot(s,u);

hold on; plot(s,v); hold off;



title ("Runge-Kutta");

h = legend ("u(s)", "v(s)");

legend (h, "location", "northwestoutside");

set (h, "fontsize", 10);

You misunderstood something in the method. The intermediate values for p,q are computed the same way as the intermediate values for u,v, and both are "Euler steps" with the last computed slopes, not separate slope computations. For the first ones that is

  u1(j) = u(i) + (1/2)*k1_u(j)*h;

  v1(j) = v(i) + (1/2)*k1_v(j)*h;

  p1(j) = p(i) + (1/2)*k1_p(j)*h;

  q1(j) = q(i) + (1/2)*k1_q(j)*h;

The computation for the k2 values then is correct, the next midpoints need to be computed correctly via "Euler steps", etc.