Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

函数功能编辑本段回目录

ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)³。解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解.

www.iLoveMatlab.cn

使用方法编辑本段回目录

[T,Y] = ode45(odefun,tspan,y0)

odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名 www.iLoveMatlab.cn

tspan  是区间 [t0 tf] 或者一系列散点[t0,t1,...,tf] book.iLoveMatlab.cn

y0     是初始值向量

www.iLoveMatlab.cn

T      返回列向量的时间点 www.iLoveMatlab.cn

Y      返回对应T的求解列向量

www.iLoveMatlab.cn

[T,Y] = ode45(odefun,tspan,y0,options) 《Simulink与信号处理》

options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 《Simulink与信号处理》

[T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options)

Matlab中文论坛

在设置了事件参数后的对应输出 www.iLoveMatlab.cn

TE      事件发生时间 book.iLoveMatlab.cn

YE      事件解决时间

Matlab中文论坛

IE      事件消失时间

www.iLoveMatlab.cn

sol =ode45(odefun,[t0 tf],y0...)  book.iLoveMatlab.cn

sol     结构体输出结果

www.iLoveMatlab.cn

应用举例编辑本段回目录

1 求解一阶常微分方程

程序:

一阶常微分方程

odefun=@(t,y) (y+3*t)/t^2;    %定义函数

《Simulink与信号处理》

tspan=[1 4];                 %求解区间

y0=-2;                       %初值 《Simulink与信号处理》

[t,y]=ode45(odefun,tspan,y0);

plot(t,y)                    %作图

title('t^2y''=y+3t,y(1)=-2,1

 legend('t^2y''=y+3t')

 xlabel('t') 《Simulink与信号处理》

 ylabel('y')

%  精确解

%  dsolve('t^2*Dy=y+3*t','y(1)=-2')

% ans =

一阶求解结果图

% (3*Ei(1) - 2*exp(1))/exp(1/t) - (3*Ei(1/t))/exp(1/t)

book.iLoveMatlab.cn

2 求解高阶常微分方程

关键是将高阶转为一阶,odefun的书写.

F(y,y',y''...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为列向量

dxdy=[y(1),y(2)....]

《Simulink与信号处理》

程序:

function Testode45

tspan=[3.9 4.0]; %求解区间

y0=[2 8];    %初值

[t,x]=ode45(@odefun,tspan,y0);

plot(t,x(:,1),'-o',t,x(:,2),'-*')

legend('y1','y2')

title('y'' ''=-t*y + e^t*y'' +3sin2t') book.iLoveMatlab.cn

 xlabel('t')

 ylabel('y')

function y=odefun(t,x)

y=zeros(2,1); % 列向量

y(1)=x(2);

y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t); 

end

end

高阶求解结果图

相关函数编辑本段回目录

ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现(自己写的，非直接调用)

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。

(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'）

Yn+1=Yn+h*f(Xn,Yn)

另外根据微分中值定理，存在0使得

Yn+1=Yn+h*f(Xn+th,Y(Xn+th))

这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。

利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：

K1＝f(Xn,Yn);

K2=f(Xn+h/2,Yn+(h/2)*K1);

K3=f(Xn+h/2,Yn+(h/2)*K2);

K4=f(Xn+h,Yn+h*K3);

Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6);

所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。

仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数：

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

调用的子函数以及其调用语句：

function dy=test_fun(x,y)

dy = zeros(3,1);%初始化列向量

dy(1) = y(2) * y(3);

dy(2) = -y(1) + y(3);

dy(3) = -0.51 * y(1) * y(2);

对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：

[T,F] = ode45(@test_fun,[0 15],[1 1 3]);

subplot(121)

plot(T,F)%Matlab自带的ode45函数效果

title('ode45函数效果')

[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数

subplot(122)

plot(T1,F1)%自编的龙格库塔函数效果

title('自编的   龙格库塔函数')

本文来源：https://www.2haoxitong.net/k/doc/8211fbd428ea81c758f57893.html