CODEWITHSUNDEEP
matlaboctavenumerical-methodsoderunge-kutta
R.FALLEN
R.FALLEN

Reputation: 21

Runge-Kutta Numerical Method Bad Approximation

I´m attempting to use the Runge-Kutta method to compare it to the lsode function. But it is performing rather poorly, every other method I used (forwards 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.sstatic.net/vJ6Yi.png

If anyone can point out a way to improve it or if I'm 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);

Upvotes: 0

Views: 179

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

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.

Upvotes: 2

Related Questions