matlab解常微分方程

1. ODE

常微分方程ordinary differential equation的缩写，此种表述方式常见于编程，如MATLAB中Simulink求解器solver已能提供了7种微分方程求解方法：ode45(Dormand-Prince)，ode23(Bogacki-Shampine)，ode113(Adams)，ode15s(stiff/NDF)，ode23s(stiff/Mod. Rosenbrock)，ode23t(mod. stiff/Trapezoidal)，ode23tb(stiff/TR-BDF2)。

微分方程、微分方程组

自标量 因变量 一元  多元  函数  映射 

一元：只有一个因变量

多元：有多个因变量

导数 偏导：谁对谁的导数，因变量对自变量的导数，默认或缺省自变量为t 、x  ？

一元方程  多元方程  多元方程组  n个方程解n个未知量

微分方程  一阶  高阶微分方程  一阶微分方程组

一阶常微分方程：Dx/dt + x = e^t

高阶常微分方程：d^2x/dt^2+dx/dt+x=e^2t

一阶微分方程组（多元）：

dy/dt+x=e^2t

dx/dt+2y-x=e^t

初始条件：dy/dt0=...  dx/dt0=...  y0=...  x0=... 

可以解出：y=f(t)=....    x=f(t)=....  两个方程解两个未知数（因变量）

一个N阶（多元）微分方程可以写成（分解成）N个一阶微分方程（即微分方程组）

如：

x.. + 2x. -x = u

令x.=x2; x=x1  则...

微分方程的精确解： r=dsolve('eqn1','eqn2',...,'cond1','cond2',...,'var').

数值解： [t,y]=solver('odefun',tspan,y0,options)

1. 求精确解

1.微分方程

r=dsolve('eqn1','eqn2',...,'cond1','cond2',...,'var').

该命令中可以用D表示微分符号，其中D2表示二阶微分，D3表示三阶微分，以此类推。

解释如下：eqni表示第i个微分方程，condi表示第i个初始条件，var表示微分方程中的自变量，默认为t。

1 2 3

>> dsolve('Dy=3*x^2','y(0)=2','x')    ans = x^3 + 2

 2.微分方程组

1 2 3 4 5 6 7 8 9 10

>> [x,y]=dsolve('Dx=y','D2y-Dy=0','x(0)=2','y(0)=1','Dy(0)=1')    x =    exp(t) + 1       y =    exp(t)

3.求解微分方程组在初始条件x (t = 0)= 1, y (t =0 )= 0 下的特解,并画出解函数的图像。

1 2 3 4 5 6 7 8 9 10 11 12

>> [x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')    x =    exp(t*(15^(1/2) - 1))*(15^(1/2) - 4)*((13*15^(1/2))/330 - exp(2*t - 15^(1/2)*t)*(15^(1/2)/165 + 1/22) + 1/22) - exp(-t*(15^(1/2) + 1))*(exp(2*t + 15^(1/2)*t)*(15^(1/2)/165 - 1/22) + (15^(1/2)*(15^(1/2) - 13))/330)*(15^(1/2) + 4)       y =    exp(-t*(15^(1/2) + 1))*(exp(2*t + 15^(1/2)*t)*(15^(1/2)/165 - 1/22) + (15^(1/2)*(15^(1/2) - 13))/330) + exp(t*(15^(1/2) - 1))*((13*15^(1/2))/330 - exp(2*t - 15^(1/2)*t)*(15^(1/2)/165 + 1/22) + 1/22)    >> ezplot(x,y)

1. ezplot与plot的区别

plot(x,y)以x为横坐标，y为纵坐标绘制曲线

plot(x,y1,x,y2,...)以x为横坐标值，以y1，y2...元素为纵坐标值绘制多条曲线

plot中x,y的表达式是已知的或者是形如y=f(x)的表达式

而ezplot是画出隐函数图形，是形如f(x,y)=0这种不能写出像y=f(x)这种函数的图形，explot无需数据准备，直接画出函数图形

1. 求近似解

ode求解器

求解器	问题类型	精度	何时使用
ode45	非刚性	中	大多数情况下，您应当首先尝试求解器 ode45。
ode23		低	对于容差较宽松的问题或在刚度适中的情况下，ode23 可能比 ode45 更加高效。
ode113		低到高	对于具有严格误差容限的问题或在 ODE 函数需要大量计算开销的情况下，ode113 可能比 ode45 更加高效。
ode15s	刚性	低到中	若 ode45 失败或效率低下并且您怀疑面临刚性问题，请尝试 ode15s。此外，当解算微分代数方程 (DAE) 时，请使用 ode15s。
ode23s		低	对于误差容限较宽松的问题，ode23s 可能比 ode15s 更加高效。它可以解算一些刚性问题，而使用 ode15s 解算这些问题的效率不高。 ode23s 会在每一步计算 Jacobian，因此通过 odeset 提供 Jacobian 有利于最大限度地提高效率和精度。 如果存在质量矩阵，则它必须为常量矩阵。
ode23t		低	对于仅仅是刚度适中的问题，并且您需要没有数值阻尼的解，请使用 ode23t。  ode23t 可解算微分代数方程 (DAE)。
ode23tb		低	与 ode23s 一样，对于误差容限较宽松的问题，ode23tb 求解器可能比 ode15s 更加高效。
ode15i	完全隐式	低	对于完全隐式问题 f(t,y,y’) = 0 和微分指数为 1 的微分代数方程 (DAE)，请使用 ode15i。

 

1. 求解微分方程初值问题的数值解，求解范围为区间 [0,0.5] 。

1. inline()通俗的来说就是用于定义函数，使用inline定义一个函数

给a，b，x赋值即可得到y

1 2 3 4 5 6

>> f=inline('a*x+b','a','b','x'); >> f(1,2,3)   ans =        5

 

1. 求常微分方程的数值解，MATLAB的命令格式为：
   [t,y]=solver('odefun',tspan,y0,options)
   其中solver选择ode45等函数名，odefun为根据待解方程或方程组编写的m文件名，tspan为自变量的区间[t0,tf]，即准备在那个区间上求解，y0表示初始值，options用于设定误差限制。命令格式为：
   options=odeset('reltol',rt,'abstol',at)
   rt输入相对误差，at输入绝对误差。

 

数值解是采用如有限元方法、数值逼近方法、插值方法等计算方法得到的解。只能利用数值计算的结果，而不能随意给出自变量并求出计算值。

当无法藉由微积分技巧求得解析解时，这时便只能利用数值分析的方式来求得其数值解了。实际情况下，常微分方程往往只能求解出其数值解。

在数学中，刚性方程是指一个微分方程，其数值分析的解只有在时间间隔很小时才会稳定，只要时间间隔略大，其解就会不稳定。

目前很难去精确地去定义哪些微分方程是刚性方程，但是大体的想法是：这个方程的解包含有快速变化的部分。

1. 解析解能求出x=f(t)；数值解则需要制定t的某个值，如只能解出x在t=某个常数时的值，不能解出关于t的表达式。
2.  

y=dsolve（'el，e2，..?，'c1，c2，..?，'v'）

其中’e1，e2，..'为微分方程或微分方程组；

’c1，c2，...’，是初始条件或边界条件；

’v'是独立变量，默认的独立变量是’t’；

y返回解析解。如果没有初始条件，则求出通解，如

果有初始条件，则求出特解。

用字符串表示常微分方程，自变量缺省时为t，导数用

D表示微分。y的2阶导数用D2y表示，依此类推。[T，Y，TE，YE，IE]=solver（'odefun'，tspan，y0，options）

其中solver为ode23、ode45、ode113、ode15s、ode23s、

ode23t、ode23tb函数；

odefun是函数句柄；

tspan微分定义区间；

y0为初值行矩阵；

T值是t序列（为列向量）；

Y值是微分方程的解Y在各点t的值（为列向量）；

TE表示事件发生时间，可缺省；

YE表示事件解决时间，可缺省；

IE表示事件消失时间，可缺省；

options是求解参数设置，可以用odeset在计算前设定误差，输出参数，事件等，可缺省。

 

1. 在微分方程的解析解法中，我们着重强调dsolve有很大的局限性，一般没法解决非线性微分方程，由此我们有了ODEs的数值解法
   其中微分方程的转化是解决所有问题的关键点，必须掌握了它，才能过使用ode**求解器
   在这里我们具体介绍下面的几种微分方程的数值解法
   1ODE解算器简介(ode**)   http://www.matlabsky.com/thread-528-1-1.html
   2微分方程转换  http://www.matlabsky.com/thread-534-1-1.html
   3刚性/非刚性问题(Stiff/Nonstiff)  http://www.matlabsky.com/thread-530-1-1.html
   4隐式微分方程(IDE)  http://www.matlabsky.com/thread-529-1-1.html
   5微分代数方程(DAE)  http://www.matlabsky.com/thread-531-1-1.html
   6延迟微分方程(DDE)  http://www.matlabsky.com/thread-532-1-1.html
   7边值问题(BVP)  http://www.matlabsky.com/thread-533-1-1.html
   在介绍上面其中常微分方程之前，我们有必要了解下Matlab中提供的求解ODEs的解算器ode**，知道这些函数的调用和参数意义。

 

好，上面我们把MATLAB中的常微分方程(ODE)的解算器讲解的差不多了，下面我们就具体开始介绍如何使用上面的知识吧！ 现实总是残酷的，要得到就必须先付出，不可能所有的ODE一拿来就可以直接使用，因此，在使用ODE解算器之前，我们需要做的第一步，也是最重要的一步，借助状态变量将微分方程组化成Matlab可接受的标准形式(一阶显示常微分方程)！ 如果ODEs由一个或多个高阶微分方程给出，则我们应先将其变换成一阶显式常微分方程组！ 下面我们以两个高阶微分方程构成的ODEs为例介绍如何将之变换成一个一阶显式常微分方程组。 step1.将微分方程的最高阶变量移到等式的左边，其他移到右边，并按阶次从低到高排列，假如说两个高阶微分方程最后能够显式的表达成如下所示： 我们说过现实总是残酷的，有时方程偏偏是隐式的，没法写成上面的样子，不用担心Matlab早就为我们想到了，这个留在后面的隐式微分方程数值求解中再详细讲解！ step2.为每一阶微分式选择状态变量，最高阶除外 从上面的变换，我们注意到，ODEs中所有因变量的最高阶次之和就是需要的状态变量的个数，最高阶的微分式(比如上面的x (m)和y(n))不需要给它一个状态变量 step3.根据上面选用的状态变量，写出所有状态变量的一阶微分的表达式 注意到，最高阶次的微分式的表达式直接就是step1中的微分方程 好，到此为止，一阶显式常微分方程组，变化顺利结束，接下来就是Matlab编程了。请记住上面的变化很重要，Matla中所有微分方程的求解都需要上面的变换。 下面通过一个实例演示ODEs的转换和求解 【解】真是万幸，该ODEs已经帮我们完成了step1，我们只需要完成step2和step3了 (1)选择一组状态变量 (2)写出所有状态变量的一阶微分表达式 (4)有了数学模型描述，则可以立即写出相应的Matlab代码了(这里我们优先选择ode45) x0=[1.2;0;0;-1.04935751];%x0(i)对应与xi的初值 options=odeset('reltol',1e-8);%该命令的另一种写法是options=odeset;options.reltol=1e-8; tic [t,y]=ode45(@appollo,[0,20],x0,options);%t是时间点，y的第i列对应xi的值，t和y的行数相同 toc plot(y(:,1),y(:,3))%绘制x1和x3，也就是x和y的图形 title('Appollo卫星运动轨迹') xlabel('X') ylabel('Y')   function dx=appollo(t,x) mu=1/82.45; mustar=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mustar)^2+x(3)^2); dx=[x(2)     2*x(4)+x(1)-mustar*(x(1)+mu)/r1^3-mu*(x(1)-mustar)/r2^3     x(4) -2*x(2)+x(3)-mustar*x(3)/r1^3-mu*x(3)/r2^3];

 

 

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

 

1. 例子