[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];
4
2.微分方程转换
好， 上面我们把 Matlab 中的常微分方程(ODE)的解算器讲解的差不多了， 下面我们就具体开始介绍如何使用上面 的知识吧！ 现实总是残酷的，要得到就必须先付出，不可能所有的 ODE 一拿来就可以直接使用，因此，在使用 ODE 解算 器之前，我们需要做的第一步，也是最重要的一步就是将微分方程组化成 Matlab 可接受的标准形式！ 如果 ODEs 由一个或多个高阶微分方程给出，则我们应先将其变换成一阶显式常微分方程组！ 下面我们以两个高阶微分方程构成的 ODEs 为例介绍如何将之变换成一个一阶显式常微分方程组。
【解】真是万幸，该 ODEs 已经帮我们完成了 step1，我们只需要完成 step2 和 step3 了
(1)选择一组状态变量
x1 = x, x2 = x ', x3 = y, x4 = y '
(2)写出所有状态变量的一阶微分表达式
6
x1 ' = x2 x2 ' = 2 x4 + x1 − µ* ( x1 + µ ) / r13 − µ ( x1 − µ * ) / r23 x3 ' = x4 x4 ' = −2 x2 + x3 − µ * x3 / r13 − µ x3 / r23
在求解 ODE 时，我们还会用到 deval()函数，deval 的作用就是通过结构体 solution 计算 t 对应 x 值，和 polyval 之类的很相似 输入参数： sol：就是上次调用 ode**函数得道的结构体解 xint：需要计算的点，可以是标量或者向量，但是必须在 tspan 范围内 该函数的好处就是如果我想知道 t=t0 时的 y 值，不需要重新使用 ode 计算，而直接使用上次计算的得到 solution 就可以
注意到，最高阶次的微分式的表达式直接就是 step1 中的微分方程 好，到此为止，一阶显式常微分方程组，变化顺利结束，接下来就是 Matlab 编程了。请记住上面的变化很重要， Matla 中所有微分方程的求解都需要上面的变换
下面通过一个实例演示 ODEs 的转换和求解
已知Apoll卫星的运动轨迹( x, y )满足下面的微分方程组 µ *( x + µ ) µ (x − µ* ) − x '' = 2 y '+ x − r13 r23 * y '' = −2 x '+ y − µ y − µ y r13 r23 其中： r1 = ( x + µ )2 + y 2 , r2 = ( x − µ * )2 + y 2 , µ * = 1 − µ , µ = 1/ 82.45 x(0) = 1.2, x '(0) = 0, y (0) = 0, y '(0) = −1.04935751
非刚性
低到高
刚性(stiff) 刚性 中等刚性 刚性
低到中 低 低 低
相关优化选项 参数 relto abstol normcontrol outputfcn outputsel refine stats nonnegative events maxstep inialstep jacobian jpattern vectorized mass mstatedependence ode45 √ ode23 √ ode113 √ ode15s √ ode23s √ ode23t √ ode23tb √
Matlab 中常微分方程数值解法讲解©
1.ODE 解算器简介 ..............................................................................................................................................................3 2.微分方程转换 ...................................................................................................................................................................5 3.刚性/非刚性问题 ..............................................................................................................................................................8 4.隐式微分方程(IDE) ........................................................................................................................................................10 5.微分代数方程(DAE).......................................................................................................................................................15 6.延迟微分方程(DDE).......................................................................................................................................................18 7.边值问题(BVP) ...............................................................................................................................................................20
然而现实总是残酷的，有时方程偏偏是隐式的，没法写成上面的样子，不用担心 Matlab 早就为我们想到了，这 个留在后面的隐式微分方程数值求解中再详细讲解！
step2.为每一阶微分式选择状态变量，最高阶除外
x1 = x, x2 = x ', x3 = x '', L, xm = x ( m −1) xm+1 = y, xm+1 = y ', xm + 2 = y '',L , xm + n = y ( n−1)
step1.将微分方程的最高阶变量移到等式的左边，其他移到右边，并按阶次从低到高排列，假如说两个高阶微
分方程最后能够显式的表达成如下所示：
x ( m ) = f (t , x, x ', x '', L , x ( m −1) , y , y ', L , y ( n−1) ) (n) ( m −1) , y , y ', L , y ( n−1) ) y = g (t , x, x ', x '', L , x
√ √ √ √
× × × ×
√ √ √ ×
√ × × ×
注：√*表示只能应用于没有质量矩阵的问题 输入参数： odefun：微分方程的句柄，必须写成 Matlab 规范格式，这个具体在后面讲解 tspab=[t0 tf]或者[t0,t1,…tf]：微分变量的范围，两者都是根据 t0 和 tf 的值自动选择步长，后只是前者返回所有计 算点的微分值，而后者只返回指定的点的微分值，一定要注意对于后者 tspan 必须严格单调，还有就是两者数据 存储时使用的内存不同(明显前者多)，其它没有任何本质的区别 y0=[y(0),y'(0),y"(0)…]：微分方程初值，依次输入所有状态变量的初值 options：微分优化参数，是一个结构体，使用 odeset 可以设置具体参数，详细内容查看帮助 输出参数： T：时间列向量 Y：二维数组，第 i 列表示第 i 个状态变量的值， 行数与 T 一致
从上面的变换，我们注意到，ODEs 中所有因变量的最高阶次之和就是需要的状态变量的个数，最高阶的微分式 (比如上面的 x-(m)和 y(n))不需要给它一个状态变量
step3.根据上面选用的状态变量，写出所有状态变量的一阶微分的表达式
5
x '1 = x2 x '2 = x3 x '3 = x4 M x 'm = f (t , x1 , x2 , x3 ,L , xm + n ) x 'm+1 = xm + 2 M x 'm+ n = g (t , x1 , x2 , x3 ,L , xm+ n )