用数值方法求解微分方程组

本楼主要讲解简单的Matlab解初值问题的方法，利用odesolver来进行求解。

刚性方程
如果在一个过程中的快变量过程与慢变量过程的变化速度相差非常大，在数学上称这种过程具有刚性(个人理解是对参数变化的敏感程度)。描述这宗过程的微分方程组称为刚性方程组。也称为病态方程组或坏条件方程组。否则称之为非刚性方程组。对于一阶常系数线性微分方程组

若其雅克比矩阵A的特征值相差十分悬殊，称此微分方程组是刚性的。

odesolver 简要介绍
odeij 问题类型 精度 适用对象
ode45 非刚性 中度 多数情况优先适用，但不能用来解决刚性问题
ode113 非刚性 低到高 不能解决刚性问题，当误差容限要求严格时，效果较ode45好
ode23 非刚性 较低 用来解决刚性问题，或允许误差范围较宽的问题。
ode15s 刚性 低到高 可用来解决刚性问题，当采用ode45失败或效果很差时，可考虑采用。
ode23s 刚性 低 可用来解决刚性问题，误差容限较宽时效果较ode15s好。
ode23tb 刚性 低 可用来解决刚性问题，误差容限较宽时效果较ode15s好。
ode23t 适度刚性 低 可用来解决刚性问题，但要求无数值衰减
ode23，ode45等库函数主要采用自适应龙格库塔方法。其中ode23系列采用二阶、三阶龙格库塔方法，ode45系列采用四阶、五阶龙格库塔方法。一般ode45积分段少，运算速度更快

格式调用：
1. [T, Y] = odeij ( @odefun,tspan, y0)
2. [T, Y] = odeij ( @odefun,tspan, y0, options )
3. [T, Y] = odeij ( @odefun,tspan, y0, options, P1, P2, … )
输入量：
odefun是待解的微分方程组中函数f( t, y )写成m函数的文件名。
tspan = [t0 tfinal]是积分区间。其中t0、tfinal分别为自变量的初值与终值。如果要求输出在指定点t=t0, t1, t2, …的微分方程数值解y=y0, y1, y2, …则输入tspan = [ t0 t1 t2 … ]；等步长用t=t0 : k : tfinal。
t的步长由系统自动选取
options是用来设置可选参数的。一般设置：RelTol相对误差(default = 1e-3)与AbsTol绝对误差(default = 1e-6)
options =odeset()…
P1, P2, …的作用是将附加参数P1、P2…传递给ode文件(如odefun( t,y,P1, P2…))和在
options中的所有函数。如果没有options参数设置，则用options=[]代替，以便正确传递函数。如果tspan或y0没有定义，则odeij函数调用ode文件，令
[tspan, y0, options] = ODEFUN ( [], [],’init’ )来获得参数，甚至调用参数列表中末尾的没有的参数。

举一个最简单的例子

精确解：
y=dsolve('Dy-f(x,y)=0','y(x0)=y0','x')
数值解：
1. 编写文件funfcn.m
function f = funfcn( x, y )
f = f( x, y )
2. 保存，并运行
options = odeset( ‘Reltol’, …, ‘AbsTol’, … )
或options=odeset; options.RelTol=…; options.AbsTol=…;
[ t, y ] = ode23( @funfcn, [a, b], y0,options )

数值解：
1. 编写文件funfcn.m
function f = funfcn( x, y )
f(1) = f( x, y ,z)
f(2) = g(x,y,z)
2. 保存，并运行
options = odeset( ‘Reltol’, …, ‘AbsTol’, … )
或options=odeset; options.RelTol=…; options.AbsTol=…;
[ t, y ] = ode23( @funfcn, [a, b], [y0,z0],options )

ode函数在控制领域中应用广泛，主要是一些微分方程难以求解，同时工程上，数值的意义远大于表达式的意义。

简要介绍一下ode函数调用自定义函数，自定义函数建立的方法。

在Command window中就可以使用
[t,x]=ode45('fun',[0,tfinal],x0)%取tfinal=100,x0=[2,3]
对于两个变量x和y，还可以绘出
相平面曲线
plot(x(:,1),x(:,2))

换一个复杂一点的

其中含有参量Param。

这个时候就需要第三个调用格式来传递参数了。
建立
function x_dot=fun(t,x,flag,alpha,beta,lambda,a,b)
x_dot=[alpha*(x(2)-x(1)-(b*x(1)+(a-b)/2*(abs(x(1)+1)-abs(x(1)-1))));
x(1)-x(2)+x(3);
-beta*x(2)-lambda*x(3)];
或
fun=inline(['[alpha*(x(2)-x(1)-(b*x(1)+(a-b)/2*(abs(x(1)+1)-abs(x(1)-1))));',...
'x(1)-x(2)+x(3);-beta*x(2)-lambda*x(3)]'],'t','x','flag','alpha','beta','lambda','a','b');
t_final=0.15;
x0=[-2.121304;-0.066170;2.881090];
alpha=15;beta=20;lambda=0.5;a=-120/7;b=-75/7
其中inline函数也可以建立自定义函数，语法类似，这里不再讲述。
输入
[t,x]=ode45('fun',[t0,tfinal],x0,[],alpha,beta,lambda,a,b)
绘出图像：
plot(t,x);axis([0,.15,-1e5,1e5])
绘出相空间曲线
plot3(x(:,1),x(:,2),x(:,3));

以上可以看出，当微分方程形式为
时，建立function函数x_dot，参数包括
t：自变量
x：因变量
flag：用于占位，可能可以省略吧，忘了
PARAM 1，PARAM 2，PARAM 3，… ：方程组中的附加参数，可以方便的在不改变m文件的情况下，调整各个参数值。
建立的列矩阵x_dot
每行为对应单个方程右端的式子，其中变量由向量替换
同样，function的m文件也可由inline函数创建，规则如例中所示。

高阶微分方程组的求解
学过差分方程的都知道，高阶微分方程可以转换为微分方程组的形式。