MATLAB求解常微分方程初值问题的三类方法

求微分方程初值问题：

我们以二阶微分方程为例来介绍这三种方法的使用，当然，第一种方法——四阶龙格-库塔法实际上就是自己手动编写程序来求解，而第二种方法——ode45就是matlab上的也是用四阶龙格-库塔法实现的更完备的方法，都是数值解法，第三种方法是微分方程的符号解法，是精确解法。

讨论一阶微分方程

初始条件为，其中y和f可以是 矢量 （即微分方程组）。

所谓的数值积分，是 解决 在初值已知的情况下，对进行近似积分的问题，通常称为微分方程的初值问题。

对上述一阶微分方程两边进行积分：

在离散的时间下：

注：如果要求解的是高阶微分方程，可降阶为一阶微分方程组。如二阶微分方程 可降阶为： 也就是把原来二阶的看成：，将一阶导数令为：

步长： 单步法：只由前一时刻的数值就可以求得后一时刻的数值，是一种自启动算法。（如欧拉法、龙格-库塔法等） 多步法：计算时要用到y的数值，这种方法为多步法，不是自启动算法。（如Adams法） 截断误差：分析数值积分的精度常用泰勒级数，设准确，则: 若只取前几项之和来近似计算，这样产生的误差称为截断误差，比如取到，就称具有k阶精度，即该方法是k阶的。

一、四阶龙格-库塔法

在看四阶之前，先了解一下二阶的龙格-库塔法。

二阶，即泰勒展开到二阶：

将

代入上式，然后假设微分方程的解可写成如下形式：

取不同的可以得到不同的二阶龙格-库塔法。

那么四阶龙格-库塔法 类  似：

具体计算时，由可以求得，由可以求得，……，直到所要求的值。

附上matlab 代码框架 ：

首先定义一阶微分方程组 ：

function func=f(t,y)f1=y(2);f2=(1-y(1)^2)*y(2)-y(1);func=[f1;f2];end

然后是龙格-库塔算法：

function [t,Y]=RKutta(Delta,Y0,tm)% Delta为步长，Y0为初始值，是一个列向量,tm是计算Y(:,1)=Y0;k=1;t=0:Delta:tm;t_size=length(t);for i=1:t_size    z1=f(t(k),Y(:,k));    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;    k=k+1;endY=Y(:,1:end-1); % 其实求得的是0到20+Delta秒的yend

设定初值和步长，进行计算：

Delta=0.001;Y0=[2;0];tm=20;[t1,y1]=RKutta(Delta,Y0,tm)

将y随时间变化规律画出来：

二、ode45

matlab提供了好几种求解器，其中最常用的就是ode45求解器。

其最简单的 使用方法 是：[t,y] = ode45(odefun,tspan,y0)

同样首先定义求解的微分方程组odefun：

function dydt = vdp1(t,y) dydt = [y(2); (1-y(1)^2)*y(2)-y(1)]; end

求解时间范围tspan=[0,20]，初值y0=[2;0]

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

求得的是y随着t变化的序列。

同样画出图像：

此外，如果积分区间不定，但知道需要在y减小到-1.5时停止积分，就需要使用ode45的事件描述：

function [position,isterminal,direction] = myEventsFcn(t,y)position = y(1)+1.5; % 该式等于0即事件发生（即y(1)=-1.5）isterminal = 1;  % 当事件多于1个时，设置当某个事件发生时停止积分（为1）direction = -1;   % 值为 -1 时仅在事件函数递减的位置查找零

ode45的 事件 设置：

options=odeset('events',@myEventsFcn);[t2,y2] = ode45(@vdp1,[0 20],[2; 0],options);

此时画出图像为：

三、dsolve

符号解也是精确解，然而大多数方程并不能求得解析解，比如本例中的这个二阶微分方程，用此方法求解会得到报错：

syms t y(t)eqn=diff(y,t,2)-(1-y^2)*diff(y,t,1)+y==0;Dy = diff(y,t);cond=[y(0)==2,Dy(0)==0];ySol(t)= dsolve(eqn,cond)

所以还是需要用数值方法去求解。

这里举一个可以求解精确解的例子：

syms t y(t)eqn=diff(y,t,2)-3*diff(y,t,1)+2*y==1;Dy = diff(y,t);cond=[y(0)==1,Dy(0)==0];ySol(t)= dsolve(eqn,cond)

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删