Runge-Kutta数值方法差逼近

Question

我尝试使用Runge-Kutta方法将其与lsode函数进行比较。但它的性能相当差，我使用的其他方法(向前和向后Euler，Heun)都比lsode做得更好，以至于它们与lsode几乎没有区别。

这是我的代码返回的https://i.stack.imgur.com/vJ6Yi.png

如果有人能指出改善它的方法，或者如果我做错了什么，我将不胜感激。

下面是我使用的Runge-Kutta方法

%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);

Answer 1

您误解了该方法中的某些内容。计算p,q的中间值的方法与计算u,v的中间值的方法相同，并且这两个值都是具有上次计算的坡度的“欧拉步长”，而不是单独的坡度计算。对于第一个，即

  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;

因此，k2值的计算是正确的，需要通过“欧拉步”正确地计算下一个中点，等等。