6. MATLAB解微分方程问题(龙格库塔法)

最新推荐文章于 2024-06-25 20:21:35 发布

洛阳八中我最棒

最新推荐文章于 2024-06-25 20:21:35 发布

阅读量4k

点赞数 4

CC 4.0 BY-SA版权

分类专栏：数学建模文章标签： python matlab java 人工智能 js

本文链接：https://blog.youkuaiyun.com/weixin_44532671/article/details/110466433

本文介绍了如何用MATLAB的龙格库塔法解决微分方程问题，包括lorenz系统和一个具体的作业实例。在作业一中，针对n=2的情况进行了求解，虽然尝试调整参数，但未能达到预期的运行结果。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

老师说系统给的ode45好多都解决不了。

1.lorenz系统

test.m

clear; 
clc;

%系统龙格库塔法
[t,h] = ode45(@test_fun,[0 40],[12 4 0]);
plot3(h(:,1),h(:,2),h(:,3));
grid on;

%自定义龙格库塔法
[t1,h1]=runge_kutta(@test_fun,[12 4 0],0.01,0,40);
figure;
plot3(h1(1,:),h1(2,:),h1(3,:),'r')
grid on;

runge_kutta.m

%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点(参数形式参考了ode45函数)
function [x,y]=runge_kutta(ufunc,y0,h,a,b)
n=floor((b-a)/h);       %步数
x(1)=a;                 %时间起点
y(:,1)=y0;              %赋初值，可以是向量，但是要注意维数
for i=1:n               %龙格库塔方法进行数值求解    
    x(i+1)=x(i)+h;    
    k1=ufunc(x(i),y(:,i));  
    k2=ufunc(x(i)+h/2,y(:,i)+h*k1/2);    
    k3=ufunc(x(i)+h/2,y(:,i)+h*k2/2);   
    k4=ufunc(x(i)+h,y(:,i)+h*k3);   
    y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;
end