方程组 Runge-Kutta 使用 matlab 的方程组的 4 阶

Question

我需要使用 Runge-Kutta 方法 4 阶来做 matlab 代码来解决方程组，但在每次尝试中我都遇到了问题并且无法解决 导数为 (d^2 y)/dx^(2) +dy/dx-2y=0 ， h=0.1 Y(0)=1 ， dy/dx (0)=-2

{clear all, close all, clc
%{
____________________TASK:______________________
Solve the system of differential equations below 
in the interval 0<x<1, with stepsize h = 0.1.
y= y1                 y(0)=0
y3= 2y1-y2            y2(0)=-2   
_______________________________________________
%}
h = 0.1;
x    = 0:h:1
N = length(x);
y   = zeros(N,1);
y3    = zeros(N,1);
g = @(x, y, y1, y2) y1;
f = @(x, y, y1, y2) 2*y1-y2;
y1(1) = 0;
y2(1) =-2;
for i = 1:(N-1)
k_1 = x(i)+y(i)
k_11=g(x(i),y,y(i))
k_2 = (x(i)+h/2)+(y(i)+0.5*h*k_1)
k_22=g((x(i)+0.5*h),y,(y(i)+0.5*h*k_11))
k_3 = (x(i)+h/2)+(y(i)+0.5*h*k_2)
k_33=g((X(i)+0.5*h),y,(y(i)+0.5*h*k_22));
k_4 = (x(i)+h)+(y(i)+h*k_33)
k_44=g((x(i)+h),y,(y(i)+k_33*h));

y3(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h  
y3(:,i)=y;

end
Answer_Matrix = [x' y3 ];}

Answer 1

你使用了函数，这不是真正必要的，但这样可能更容易更清楚地看到公式。但是，在函数中，列出函数中使用的参数。这不是必需的，并且会产生不必要的开销。

在初始条件下，您应该使用y和y3，因为那是您在循环中使用的。同样在第一个条件下，你犯了一个错字。

在循环中，您忘记调用函数f，并更新y向量。

在代码中进行这些更改会导致以下结果：

h = 0.1;
x = 0:h:1;
N = length(x);
y = zeros(N,1);
y3 = zeros(N,1);
g = @(y2) y2;
f = @(y1, y2) 2*y1-y2;
y(1) = 1;
y3(1) = -2;
for i = 1:(N-1)
k_1 = f(y(i), y3(i));
k_11 = g(y3(i));
k_2 = f(y(i)+0.5*h*k_1, y3(i) +0.5*h*k_11);
k_22 = g((y3(i)+0.5*h*k_11));
k_3 = f(y(i)+0.5*h*k_2, y3(i) +0.5*h*k_22);
k_33 = g((y3(i)+0.5*h*k_22));
k_4 = f(y(i)+h*k_3, y3(i) +h*k_33);
k_44 = g((y3(i)+h*k_33));
y3(i+1) = y3(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h  ;
y(i+1) = y(i) + (1/6)*(k_11+2*k_22+2*k_33+k_44)*h  ;
end
Answer_Matrix = [x' y];
% solution of DE is exp(-2x) and is plotted as reference
plot(x,y,x,exp(-2*x))

---

如前所述，您也可以在不使用函数的情况下解决此问题：

h = .1;
x = 0:h:1;
N = length(x);
% allocate memory
y = zeros(N,1);
z = zeros(N,1);
% starting values
y(1) = 1;
z(1) = -2;
for i=1:N-1
ky1 = z(i);
kz1 = -z(i) + 2*y(i);
ky2 = z(i) + h/2*kz1;
kz2 = -z(i) - h/2*kz1 + 2*y(i) + 2*h/2*ky1;
ky3 = z(i) + h/2*kz2;
kz3 = -z(i) - h/2*kz2 + 2*y(i) + 2*h/2*ky2;
ky4 = z(i) + h*kz3;
kz4 = -z(i) - h*kz3 + 2*y(i) + 2*h*ky3;
y(i+1) = y(i) + h/6*(ky1 + 2*ky2 + 2*ky3 + ky4);
z(i+1) = z(i) + h/6*(kz1 + 2*kz2 + 2*kz3 + kz4);
end
% exp(-2*x) is solution of DE and is plotted as reference 
plot(x,y,x,exp(-2*x))