在系统中使用 Runge-Kutta 积分方法

Question

h=0.005;                                            
x = 0:h:40;                                         
y = zeros(1,length(x)); 
y(1) = 0;                                         
F_xy = ;                    

for i=1:(length(x)-1)                              
    k_1 = F_xy(x(i),y(i));
    k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
    k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));
    k_4 = F_xy((x(i)+h),(y(i)+k_3*h));

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

我有下面的代码，我觉得是对的。我知道 F_xy 上缺少部分，因为这是我的后续问题。

我有 dx/dt = = -x(2 - y) 其中 t_0 = 0, x(t_0) = 1

并且 dy/dt = y(1 − 2x) 其中 t_0 = 0, y(t_0) = 2。

我的问题是我不知道如何将这些方程式输入代码。感谢所有帮助

Answer 1

您以不一致的方式同时使用t 和x 作为自变量。从实际的微分方程来看，自变量是t，而二维系统的因变量是x和y。它们可以组合成一个状态向量u=[x,y] 然后将系统编码为接近您编写的一种方法是

h=0.005;                                            
t = 0:h:40;   
u0 = [1, 2]                                       
u = [ u0 ]                                         
function udot = F(t,u)
     x = u(1); y = u(2);
     udot = [ -x*(2 - y),  y*(1 - 2*x) ]
end                 

for i=1:(length(t)-1)                              
    k_1 = F(t(i)      , u(i,:)        );
    k_2 = F(t(i)+0.5*h, u(i,:)+0.5*h*k_1);
    k_3 = F(t(i)+0.5*h, u(i,:)+0.5*h*k_2);
    k_4 = F(t(i)+    h, u(i,:)+    h*k_3);

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

有解决方案输出

Answer 2

F_xy是你的导函数吗？

如果是这样，只需将其编写为辅助函数或函数句柄即可。例如，

F_xy=@(x,y)[-x*(2-y);y*(1-2*x)];

另请注意，您的k_1, k_2, k_3, k_4, y(i) 都是二维的。您需要重新调整 y 的大小并相应地在迭代步骤中重写索引。