ODE45函数的使用——翻译

在Matlab 中使用ode45简介
Matlab 中常微分方程常用的函数是ODE45，这个函数能够利用——龙哥库塔法-—有效求解带时间变量步长的计算。

Ode45用于求解如下的一般问题:
)(()00,,x t x x t f dt
dx ==(1） 其中，时间t 是独立变量，x 为时间相关矢量，)(x t f ,是时间t 和x 的函数.当(1）右边的)(x t f ,是固定的,且给定x 的初始值，那么问题的解是唯一的。

在ME175中，解法是不完整的，但是只要你解决了问题，就可以获得ODE 代表的系统运动趋势.这有利于得到一个直观的印象,看起来很复杂的常微分方程，代表的质点运动轨迹确实简单明了的.以下简要解释如何得到运动轨迹:
第一步:
对给定的ODE 方程进行降阶处理，得到一系列一阶方程
这就是你要做的第一步，在一张草稿纸上处理.例如，给定ODE 方程如下：
1,3,5002-===-+⋅
⋅⋅⋅y y y e y y m y （2）
对本问题，矢量x 有两个组成分量：y 和⋅y ,或 ()()⋅==y
x y
x 21（3）
且 ()()()()()()()()()()()
211251221x e x m dt x d x dt x d x +-==（4） 其中，用（3）中的式子代表了y ，⋅y ,⋅
⋅y ,于是把（2）改写为（4）.
如果求解的问题有更多阶数更多变量呢?例如,我们除了有上面的方程（2），同时还有以下的方程:
().1,0,sin 002233===-+⋅z z t z dt z d dt z d （5） 那么，我们可以通过构造更大的矢量x 同时求解y ，z ：
()()()⋅⋅⋅
===z
x z x z
x 543 （6）
然后
⎥⎦
⎤⎢⎣⎡=⋅⋅⋅⋅z z z y y x ,,,, (7） 以及 [⎥⎦
⎤==⋅⋅=⋅==⋅==000000,,,,t t t t t t z z z y y x （8） 其中,y 变量和z 变量的放置位置对求解不造成影响。

实际上，任意次序都是有效的，例如
⎥⎦⎤⎢⎣⎡=⋅⋅⋅⋅z z y z y x ,,,, 和⎥⎥⎦⎤⎢⎣⎡=⋅
⋅⋅⋅⋅⋅y y z z z x ,,,, 但是重要的是，在整个计算过程中，你使用的顺序都必须和一阶ODE 方程中定义的变量顺序相同。

之后，如果你使用的是(7)中给定的的式子,那么系统的一阶ODE 方程，由以下方程组组成。

（10）
而涉及的表征变量⎥⎦
⎤⎢⎣⎡=⋅⋅⋅⋅z z y z y x ,,,,结果如下： （11）
基本上，可以处理任意数量的高阶ODE 方程.重要的是把它们处理成多个一阶的ODE 方程，并且确保记住被求解的矢量X 中,不同变量所分配的顺序。

第二步 编写代码
既然你已经有所求解问题的一阶格式,在你编程的主要代码中，将会用到以下的命令
[]()options xinit tspan fname ode x t ,,,@45,=
·fname 是函数的M 文件名用于求解方程（1）右边代数式的值。

这个函数将被输入一阶ODE 系统中,并且被积分（见（10），（11)）。

后面,将会更详细的解释。

注：当然关于ODE45如何积分给定的方程有细微的差别,但是对于简单的问题，不分先后次序的积分,是可以接受的。

·tspan 是矢量定义了积分的起始点和终点,同时也定义了时间步长.例如,我们需要积分t=0到t=10,希望步数是100步，那么tspan=［0:0。

1：10]或者tspan=linspace(0，10，100). ·xinit 是初始条件矢量.确保初始值的顺序和给定的x 中变量和它倒数的顺序是一致的。

同时注意如果x 有5个变量,那么同时要输入5个初始值。

·option 这个在matlab 的帮助文件中有很好的说明。

对于大部分的问题,使用默认值就可以满足计算要求.
·t 是独立变量，计算数组x 在时间点t 的数值.这个矢量不必等于tspan ，ODE45自动调节步数以取得最大的效率和精确度.(在快速变化部分采用小步长，在变化缓慢部分采用大步长).
·x 相关内容如下。

X 是数组或矩阵，大小为length （t ）＊length （xinit ）.每一列x 代表不
同的因变量。

例如,⎥⎦
⎤⎢⎣⎡=⋅⋅⋅⋅z z z y y x ,,,,，为简单假定t=0,1，2.。

.。

，10,将会计算函数在11个点的值.
(12)
如果⋅z 是x 的第四个变量,那么()4,1x 得到了⋅z 在t=0时候的值，()4,7x 得到了⋅
z 在t=6时候的值，()4,11x 得到了t=10时候的值。

简而言之，
·()k x :,代表x 的第k 个变量,k=1与变量y 相关,k=2与变量⋅y 相关。

·():,j x 计算所有变量在某一时间点j 的数值 注：在产生hokey pokey 舞蹈前，史前儿童围坐在篝火前齐唱：
You put your left foot in
You put your left foot out ’
You put your left foot in
And you shake it all about
当你使用matlab 函数ODE45及时完成作业时，x 就是要做得全部内容。

不幸的是由于缺乏matlab 软件，使得这本书过时了.
·()()时间，变量值x k j x ≡,
命令[]()options xinit tspan fname ode x t ,,,@45,=的作用是重新定义变量。

、，然后，如果使用变量顺序为
⎥⎦
⎤⎢⎣⎡=⋅⋅⋅⋅z z z y y x ,,,,，应该这么写程序:
[]()options xinit tspan fname ode x t ,,,@45,=
()
()
()()
()
5:,4:,3:,2:,1:,x zdotdot x zdot x z x ydot x y =====
当然，也不应该认为定义y ，ydot 麻烦.直接表达为x 的形式（例如,使用()1:,x 代表y ），清晰的定义方式有利于后面的调试。

以下,你将以(或是被要求）绘图的形式描述感兴趣的轨迹：质点随时间运动轨迹，在平面中表示角度和径向关系等。

绘图和绘制子图的命令在matlab 帮助文件中有清晰的说明，这里不再详细说明。

记住如果你想在一张图中放多个图，应该使用子图的概念,当然在一张图片中画多个图，不是一个好主意。

别忘记给图加标签：包括标题，x 轴，y 轴的含义，如果多条曲线应该分别标明。

最后，请注意在[]()options xinit tspan fname ode x t ,,,@45,=中包含的仅仅是变量而已，依据自己的喜好使用字母，T 替换t,x0替换xinit 都是可以的.只有记住使用新变量名，之后的每个引用都用一样的名称。

另一个普遍的错误在于，同一变量的重复定义。

例如，定义()5:,x x =，如果足够幸运的话,会有错误警告;不幸的话，这种错误很难发现,要花数小时时间检查您的程序以解决问题。

此外,fname 是什么呢？
回忆下，我们还没有告诉matlab 程序应该对什么函数进行积分，是吧?这就是为什么需要fname 文件，fname 文件含有所有之前在稿纸上重写的ODE 一阶函数。

你可以对这个文件起任意的名称，只要与[]()options xinit tspan fname ode x t ,,,@45,=中使用的fname 一致。

例如，你对fname 取名superman 那么
[]()options xinit tspan erman ode x t ,,,sup @45,=是对的,而
[]()options xinit tspan batman ode x t ,,,@45,=就不正确了。

更进一步说，函数不必须像在原始代码中写的在同一个m 文件中，一些人喜欢在程序末尾书写子程序，特别是代码不长,比较简单的时候。

例如,你的代码名称ME175example 文件，那么m 文件将如下：
Function dxdt=ME175example （t,x)
％这里 t,x 和 dxdt are 只是变量而已。

你可以起任意名称
％只要 t 是独立变量，而X 是因变量
%dxdt 是推到的一阶因变量
% 定义常数
m= 1；
%定义变量使之清晰易懂
%Recall that x = [y ， ydot, z, zdot, zdotdot]y=x （1）；
ydot=x(2）；
z=x(3)；
zdot=x （4）；
zdotdot=x(5）;
%注意x 仅仅是1列五行的数组
％［t ，x ］ = ode45（＠fname ， tspan, xinit, options)
％数组dxdt 与x 的大小相同
dxdt = zeros （size(x))；
dxdt （1) = ydot ；
dxdt （2） = 1/m （5 x （2)exp （y ） +y2) ; %This is ydotdot
dxdt （3） = zdot;
dxdt(4） = zdotdot;
dxdt （5） = t-zdotdot+sin(z ）; %This is zdotdotdot
%Note that the input arguments must be t and x （in that order ） even in the case where t is not explicitly used in the function 。

基本模板
以下是基本模板,当你想对一个高阶常微分方程进行积分时，把它复制黏贴到Matlab 中
Function 任何你想要的名字
％定义起始时间tstat ，终止时间tend,时间步数n
Tstart=？；
Tend=？；
N=？；
Tspan=linspace （Tstart ， Tend ， N ）；
%定义初始值,确保正确的顺序
Xinit=[。

.； .。

.; .。

.； 。

.。

；。

.。

]；
%获得矢量x.把option 设置为默认值即可
[]()x init tspan @45,，积分函数，ode t x =
％定义输出变量
()
()
2:,1:,x ydot x y ==
... ％所需要画图的函数 ()()⎪⎭
⎫ ⎝⎛⎪⎭⎫ ⎝⎛⋅⋅⋅⋅⋅
⋅y y y y t y t y ,,, Subplot （？, ?, ?）
Plot(需要的图像）
%plot3()画3D 图
Title （’ ’)
Xlabel （' ')
Ylabel(’ ’）
Zlabel （’ ’）
——-—-—--—————-—-——-—--—---—--——--—-
%积分函数，可以作为独立m 文件，或在本程序底部，语法如下
Function dxdt=积分函数(t ，x ）
%定义积分中使用的常数 ％定义新变量；()
().
21etc x ydot x y ==
％也可以不定义，这样的好处是，便于其他人阅读你的程序
％写下你得到的一阶ODE 方程
dxdt=zeros(size （x)）
dxdt （1）=?；
dxdt （2)=？;
dxdt （3）=?;
Etc 。

小结:
如果所有的步骤都对了，方程也是正确的。

你将发现你的到的关系图是很优美的,不用你艰难的推到。

干的好！
最后，同样地问题又多种不同的解法。

多用matlab 尝试，你将有所发现，比起其他算法有
些算法非常有效率。

我使用Matlab三年了,仍然可以发现很多新的有效的算法。

总之，乐在其中。

Despite what you think，MA TLABis not out to get you （it has it's hands full giving me a hard time）
感谢Nur Adila Faruk Senan 伯克利加利福尼亚大学机械工程系
您写的文档，对此帮助很大。