matlab龙格库塔法调用,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