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;
end
Y=Y(:,1:end-1); % 其实求得的是0到20+Delta秒的y
end

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

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)