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）

且

()()

()()()

()()()()()()

211251

221x 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一阶函数。你可以对这个文件起

任意的名称，只要与[]()

,

@

t,

,=中使用的fname一致。

45

,

tspan

xinit

fname

ode

options

x

例如，你对fname取名superman那么

[]()

sup

,

@

45

,

,=是对的，而

options

tspan

xinit

erman

ode

x

t,

[]()

,

t,

,

@

45

,=就不正确了。

x

options

ode

xinit

batman

tspan

更进一步说，函数不必须像在原始代码中写的在同一个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 伯克利加利福尼亚大学机械工程系

您写的文档，对此帮助很大。