龙格-库塔法MATLAB
1. matlab 新建.m文件，编写龙格-库塔法求解
函
关于工期滞后的函关于工程严重滞后的函关于工程进度滞后的回复函关于征求同志党风廉政意见的函关于征求廉洁自律情况的复函
数
function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点(参数形式参考了ode45函数) n=floor((b-a)/h); %求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值，可以是向量，但是要注意维数
for ii=1:n
x(ii+1)=x(ii)+h;
k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); k4=ufunc(x(ii)+h,y(:,ii)+h*k3); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %按照龙格库塔方法进行数值求解
end
2.另外再新建一个.,m文件，定义要求解的常微分方程函数
function dx=fun1(t,x)
dx =zeros(2,1);%初始化列向量
dx(1) =0.08574*x(2)-1.8874*x(1)-20.17; dx(2) =1.8874*x(1)-0.08574*x(2)+115.16;
3,再新建一个.m文件，利用龙格-库塔方法求解常微分方程
[T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20); %求解步骤2定义的fun1常微分方程，@fun1是调用其函数指针，从0到20，间隔为1
subplot(122)
plot(T1,F1)%自编的龙格库塔函数效果
title('自编的龙格库塔函数')
grid
运行步骤3文件即可得到结果，F1为估测值
或者可以调用matlab自带函数ode45求解
方法如下:
[T,F] = ode45(@fun1,0:1:20,[17.64 37800 ]);
subplot(121)
plot(T,F)%Matlab自带的ode45函数效果
title('ode45函数效果')
grid
浙公网安备 33021202002483