利用4－阶龙格库塔法数值求解

dx/dt=x(1–2y)
dy/dt=–y(2-x)
x(0)=1000,y(0)=500
其中t:时间、x:兔子数量、y:狐狸数量
绘制解的图像并做出解释。

四阶龙格库塔方法求解一次常微分方程组_哆啦A梦PLUS的博客-CSDN博客_四阶龙格库塔法 四阶龙格库塔方法求解一次常微分方程组一、写在前面二、四阶龙格库塔方法三、使用四阶龙格库塔方法求解一次常微分方程组一、写在前面龙格库塔方法是数值求解常微分非线性方程的有利工具，计算精度较高，通过缩短步进距离和增加阶数可以进一步控制误差范围。工程上较为常用的是四阶龙格库塔算法（R-K4)，在计算收敛的情况下往往可以得到比较好的结果。二、四阶龙格库塔方法这里简单介绍一下算法的具体实现过程，不做详细的推导。其求解的问题是形如方程：y˙=f(y,t),其中t∈[t0,t1]初值y(t0)=c0\dot{y https://blog.csdn.net/HelloWorldTM/article/details/112915551?ops_request_misc=&request_id=&biz_id=102&utm_term=%E5%88%A9%E7%94%A84%EF%BC%8D%E9%98%B6%E9%BE%99%E6%A0%BC%E5%BA%93%E5%A1%94%E6%B3%95%E6%95%B0%E5%80%BC%E6%B1%82%E8%A7%A3&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduweb~default-0-112915551.nonecase&spm=1018.2226.3001.4187

例如要求解

img


代码为

clear;
clc;
close all;

h=1e-5;      %步进长度
t=0:h:1;   %生成自变量t的向量

%%创建计算结果x，y，z的数组
N=length(t);
x=ones(1,N);
y=ones(1,N);
z=ones(1,N);

%%四阶龙格库塔迭代
for i=2:N
    t_n=t(i-1);
    x_n=x(i-1);
    y_n=y(i-1);
    z_n=z(i-1);
    
    kx1=y_n+3*z_n+sin(5*t_n);
    ky1=x_n+cos(t_n);
    kz1=x_n+z_n-3*cos(3*t_n)*sin(4*t_n);
    
    kx2=(y_n+ky1*h/2)+3*(z_n+kz1*h/2)+sin(5*(t_n+h/2));
    ky2=(x_n+kx1*h/2)+cos(t_n+h/2);
    kz2=(x_n+kx1*h/2)+(z_n+kz1*h/2)-3*cos(3*(t_n+h/2))*sin(4*(t_n+h/2));
    
    kx3=(y_n+ky2*h/2)+3*(z_n+kz2*h/2)+sin(5*(t_n+h/2));
    ky3=(x_n+kx2*h/2)+cos(t_n+h/2);
    kz3=(x_n+kx2*h/2)+(z_n+kz2*h/2)-3*cos(3*(t_n+h/2))*sin(4*(t_n+h/2));
    
    kx4=(y_n+ky3*h)+3*(z_n+kz3*h)+sin(5*(t_n+h));
    ky4=(x_n+kx3*h)+cos(t_n+h);
    kz4=(x_n+kx3*h)+(z_n+kz3*h)-3*cos(3*(t_n+h))*sin(4*(t_n+h));
    
    x(i)=x_n+h/6*(kx1+2*kx2+2*kx3+kx4);
    y(i)=y_n+h/6*(ky1+2*ky2+2*ky3+ky4);
    z(i)=z_n+h/6*(kz1+2*kz2+2*kz3+kz4);  
end
%%画图
figure();
hold on;
plot(t,x,'r');
plot(t,y,'g');
plot(t,z,'b');
legend('x','y','z');
xlabel('t');
title('xyz函数图象');
hold off;