matlab中的微分方程

zhdgzhdg

Matlab中微分方程的问题
1． Matlab能够处理什么样的微分方程？
2． 可以从什么地方获得更多的指导与附加信息？
常问问题
3． 对ODE求解器的语法存在有些什么变化？
4． 如何减小ODE的阶次？
5． 如何解决时变（依赖时间的）ODEs？
6． 如何采用定时间步长（？Fixed Time Step）？
7． 如何利用随机微分方程？
例子
8． 方程系统
9． 带边界值问题（Boundary Value Problem，BVP）：管道流（？Channel Flow）
刚度（？Stiffness）
10． 什么是Stiffness？
11． 暗示（间接？Implicit）与明示（直接？Explicit）方法
12． 实例
设置
13． 如何在解微分方程的时候改变设置？
14． 哪些设置参数可以更改？
15． 如何将options当作函数利用？
微分--数方程与他们的索引（？Index）
16. Matlab中如何解决微分-代数方程系统？

第1节 Matlab能够处理什么样的微分方程？
Matlab提供了解决包括解微分方程在内的各种类型问题的函数：
1． 常规微分方程（ODEs）的初始值问题
初值问题是用MATLAB ODE求解器解决的最普遍的问题。初始值问题最典型的是对非刚性度（？nonstiff）问题应用ODE45，对刚性度（？stiff）问题采用ODE15S。（对于stiffness的解释，请参照“什么是Stiffness”一节。）
2． 微分-代数方程（DAEs）的初值问题
在那些守恒定律规定一些变量之间满足常数关系领域经常遇到这类问题。Matlab 可以用ODE15S 或者 ODE23T解决索引（index）为1的DAEs。（对于索引的解释，请参阅“DAEs与他们的索引”一章。）
3． 边界值问题（BVPs）
这种通常要求微分方程在两边都具有特殊的条件组成。尽管他们通常不象IVPs那样经常遇到，但是他们也是工程应用中比较常见的问题。可以利用函数BVP4C来解决这类问题。
4． 时延微分方程（DDEs）
这类微分方程包含了独立变量的延迟。他们在生物与化学模型这类大量的应用中遇到，可以通过DDE23来解决这类问题。
5． 偏微分方程（PDEs）
采用PDEPE可以解决一维时空的抛物面与椭圆方程的初值、边界值的问题。而那些对更加多的一般的偏微分方程感兴趣的可以利用PDE工具箱。
更多的matlab的综合应用技术的信息请参阅Solution8314。
更多的有关matlab采用的各种求解器的算法的信息请查看下面的URLs：
● ODE 函数
● BVP 函数
● DDE 函数
● PDE 函数

第2节 可以从什么地方获得更多的指导与附加信息？
可以从MATLAB Center、网站的新闻组、文件交换点可以获得一系列资料，可以进一步解释MATLAB解决各种方程（ODE，DAE，BVP，DDE）的求解器的算法和使用。你可以下载各种方程的文章与手册，他们通常带有大量的实例。
可以从 matlab自带的帮助文件的 Mathematics|Differential Equations下找到使用指导。
Cleve Moler的《Numerical Computing with MATLAB》的第七章详细讨论了OEDs的解法，并附带有大量的实例与简单的问题练习。

第3节 对ODE求解器的语法存在有些什么变化？
在MATLAB6.5（R13）中应用ODE求解器求解的首选语法是：
[t，y]=odesolver（odefun，tspan，y0，options，parameter1，parameter2，…，parameterN）；
odesolver 是你采用的求解器，例如ODE45或者ODE15S。odefun是微分方程的定义函数，所以odefun定义独立参数（典型的是时间t）的导数y‘ 以及y和其他的参数。在MATLAB6.5（R13）中，推荐使用函数句柄作为odefun。
例如，ode45（@xdot，tspan，y0），而不是用 ode45（'xdot'，tspan，y0）。
请看采用函数句柄的好处的文档：
采用函数句柄传递你定义MATLAB求解器计算的量、例如大规模矩阵或者Jacobian模式的函数。
如果你喜好采用字符串儿传递你的函数，matlab求解器将回溯匹配。
在老的matlab版本里，通过传递标志来规定求解器的状态和恰当的计算。在MATALB6.0以及其后的版本中，这就没有必要了，可以从matlab自带的文档中发现这个差别。
如果里采用的matlab的ODE求解器的老的语法，你可以看看我们FTP站点上的各种求解器的老的实例： ftp://ftp.mathworks.com/pub/doc/papers/
前面的站点包含了BVP，DAE与DDE这三个方向的采用老的语法的实例。你可以在下面的站点中找到应用ODE45与ODE23的实例：
ftp://ftp.mathworks.com/pub.mathworks/toolbox/matlab/funfun
你可以在MATLAB Center的文件交换站点查看这些例子的更新版本。

第4节 如何减小ODE的阶次？
求解一阶ODE的代码是很直接的。然而，二阶或者三阶的ODE不能够直接应用求解。你必须先将高阶的ODE改写成一阶的ODEs系统，使得它可以采用MATLAB ODE求解器。
这是一个如何将二阶微分方程改写成两个一阶微分方程以便利用MATLAB的诸如ODE45等求解器求解的例子。下面的方程组包含了一个一阶与一个二阶微分方程：
x'= - y*exp（-t/5）+y' * exp（-t/5）+1； （1）
y''= -2*sin（t）； （2）
第一步是引入一个新的变量，使得它等于具有二阶导数的自由变量的一阶导数：
z=y' （3）
对上式两边求导如下：
z' = y'' ； （4）
将（4）式带入（2）式得到如下方程：
z'= -2*sin（t） （5）
联立（1），（3）与（5）得到三个一阶微分方程：
x'= - y*exp（-t/5）+y' * exp（-t/5）+1; (1)
z=y'; （3）
z'= -2*sin（t） （5）
既然 z＝y' ，用z代替等式（1）中的y' 。而且，因为MATLAB要求所有的导数项在左边，改写等式（3）。得到如下的方程组：
x'= - y*exp（-t/5）+z* exp（-t/5）+1; (1a)
y'=z ; (6a)
z'= -2*sin（t）； （5a）
为了利用ODE45或者是MATLAB的其他的ODE求解器求解上面的方程组，需要建立一个包含这些微分方程的函数。这个函数需要两个输入：状态量与时间，返回状态的微分，建立命名为odetest.m的函数如下：

1. function xprime=odetest(t, x)
   
2. % 既然状态量以单个向量的形式输入，我们令：
   
3. ％ x（1）＝x；
   
4. ％ x（2）＝y；
   
5. ％ x（3）＝z；
   
6. 
   
7. xprime（1）＝－x（2）* exp（－t/5）＋x（3）*exp(-t/5)+1；
   
8. ％ x'= - y*exp（-t/5）+z* exp（-t/5）+1;
   
9. 
   
10. xprime（2）＝－x（3）；
    
11. ％ y'=z
    
12. 
    
13. xprime（3）＝－2×sin（t）；
    
14. ％ z'= -2*sin（t）
    
15. 
    
16. xprime＝xprime（：）；
    
17. ％ 这是为了确保返回的是个列向量

复制代码

采用ODE23或者另外的MATLAB ODE求解器求解方程系统，定义起始和停止时间以及初识的状态向量。例如：

1. t0 = 5 ; % 起始时间
   
2. tf = 20 ; ％ 停止时间
   
3. x0 = [1 –1 3] ; ％ 初识条件
   
4. [t , s] = ode23 ( @odetest, [t0 ,tf ], x0) ;
   
5. x = s (: , 1 );
   
6. y = s (: , 2 );
   
7. z = s (: , 3 );

复制代码

求解结果作图如下：

第5节 如何解决时变（Time-Dependent）ODEs？

下面是一个带有一个时变项的常规微分方程利用MATLAB ODE求解器求解的例子。时变项可以通过一个带有已知采样时间的数据集或者是一个简单的函数定义。如果时变项通过数据集定义，则这个数据集和它的采样时间可以作为ODE 求解器的附加参数传递给函数。如果时变项是通过函数定义，则这个函数在导数函数需要的时候被调用。
本例中用到的微分方程是带有一个正弦驱动项的阻尼波动（Damped Wave）方程。
y’’（t）- beta * y’（t）+ omega^2 * y（t）= A* sin ( w0 * t – theta )
MATLAB要求微分方程能够表示成如下的一阶微分方程形式：
y’（t）= B * y（t）+ f（t）
其中，y 是一个状态向量，B是一个矩阵。应用前面一节讲的技术，将这个微分方程的MATLAB里面的定义成如下：
xdot（2）＝ beta * x（2）－ omega^2 * x (1) + A * sin（w0 * t － theta）；
xdot（1）＝ x（2）；
其中，xdot＝dx /dt，x（1）＝ y，x（2）＝dy /dt。
在本例中，beta，omega，A，w0与theta需要进行定义。他们作为附件参数传递给MATLAB ODE 求解器。

例1：时变项是一个函数
建立如下的微分方程函数：

1. ％ FUN1.M： 时变微分方程例子
   
2. function xdot = fun1 ( t , x , beta , omega , A , w0 , theta )
   
3. % 时变项是 A * sin ( w0 * t – theta )
   
4. xdot ( 2 ) = beta * x (2 ) - omega ^2 * x (1) + A * sin ( w0 * t – theta ) ; %原tech－note％% 写错了, 害的我居然检查了很久才发现，符号写错了^_^
   
5. xdot ( 1 ) = x ( 2 ) ;
   
6. xdot = xdot ( : );
   
7. % 使得xdot是一列
   
8. ％ fun1.m 结束

复制代码

在MATLAB中采用如下代码调用这个函数：

1. beta = .1 ;
   
2. omega = 2 ;
   
3. A = .1 ;
   
4. w0 = 1.3 ;
   
5. theta = pi /4 ;
   
6. X0 = [ 0 1]’;
   
7. t0 = 0 ;
   
8. tf = 20 ;
   
9. options = [];
   
10. [ t , y] = ode23 ( @fun1, [t0 , tf ], X0, options, beta, omega, A, w0, theta );
    
11. plot ( t , y );

复制代码

例2：时变项是一个数据集
如下的例子需要用到interp1 命令。 建立如下的微分方程函数：

％ FUN2.M 带有数据集合的时变微分方程的例子

1. function xdot = fun2 ( t, x, beta, omega, T, P )
   
2. pt = interp1 ( T, P, t );
   
3. xdot (2 ) = beta * x (2 ) – omega^2 * x (1) + pt ;
   
4. xdot (1 ) = x (2 );
   
5. xdot = xdot ( : ); % 使得xdot是一列
   
6. ％ fun2.m 结束

复制代码

在MATLAB中采用如下的形式调用这个函数：

1. beta = .1 ;
   
2. omega = 2 ;
   
3. A = .1 ;
   
4. w0 = 1.3 ;
   
5. theta = pi /4 ;
   
6. X0 = [ 0 1]’;
   
7. t0 = 0 ;
   
8. tf = 20 ;
   
9. T = t0 – eps : .1 : tf + theta +eps ;
   
10. P = A.* sin(w0 .*T –theta );
    
11. [ t, y ] = ode23 ( @ fun2 , [ t0, tf ], X0, [ ], theta, omega, T, P );
    
12. plot ( t, y (:,1));
    
13. hold on;
    
14. plot(t,y(:,2),'g');

复制代码

因为对interp1的调用，可能使得第二个例子比第一个例子运行时间长。

第6节 如何采用定时步长（Fixed Time Step）？

MATLAB配备的常规微分方程求解器函数采用了各种方法。ODE23是基于龙格－库塔（Runge-Kutta）（2,3）积分方法，ODE45是基于龙格-库塔（4,5）积分方法。ODE113是变阶Adams-Bashforth-Mouulton PESE求解器。各种求解器和他们采用的方法详细列表请参阅MATLAB在线文档。
MATLAB通过采取迈一步，估计在这步的误差，检查其值是大于还是小于容差，然后相应地调整步长。这些积分方法是不利于采用定步长的。采用定步长算法，在当你的信号频率大于求解器的频率的时候，你就可能丢失掉一些点，因而是危险的。采用变步长算法可以确保在低频的时候采用大的步长，而在高频的时候采用小的步长。MATLAB中的ODE求解器优化了变步长算法，采用变步长的时候能够运行更快，而且显而易见的是得到的结果是更加精确的。现在在MATLAB Central站点，一些定时步长的函数可以直接利用。这些求解器有：
ODE1 一阶Euler 方法
ODE2 二阶Euler方法
ODE3 三阶龙格－库塔法
ODE4 四阶龙格－库塔法
ODE5 五阶龙格－库塔法
这些求解器可以采用如下的与否使用：
y = ODE4( odefun，tspan，y0 )；
积分在tspan所规定是时间间隔值到的时候一步一步进行处理。时间值必须是按升序或者降序的方式排列。注意步长（tspan的连续元素间的距离）并不要求必须是均匀的。如果步长是均匀的，你或许可以采用 LINSPACE创建。例如：
tspan = linspace(t0, tf, nsteps); %t0 = 0; tf = 10, nsteps=100;

第7节 如何利用随机微分方程（SDEs）？
随机微分方程是指带有随机元素的微分方程，一个典型的随机微分方程可以写成如下形式：
dX = lambda * X dt + mu * X dW;

其中，X是我们感兴趣的量，t是时间，W是一个随机变量或随机过程，lambda与mu是问题的常数参数。
一个求解SDEs的详细介绍可以在这篇文章中找到： Desmond J Higham,“An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations”，SIAM review，Vol.43,No.3。文章中在举例说明文章的观点的时候，给出了许多MATLAB中的例子。
你可以从上面的文章中得到例子列表，同时，有关在经金融领域里面的例子，可以在下面的URL中找到：
http://www.maths.strath.ac.uk/~aas96106/algfiles.html

第8节 方程系统

I． 一个常规微分方程系统包括一些微分方程，这些方程依赖于其他的方程，例如：

这是个相当简单的例子，它可能通过解析解法或者数值解法求解。通常情况下，解析方法对很多系统都是不可用的。对于线性系统，将这些方程以矩阵的形式改写，如下所示：

具体内容请参阅Gilber Strang的“Introduction to Linear Algebra and Its Application"一书获得详细的信息。
本技术手册着重在于常规微分方程的数值解法。为了用数值解法解决上面方程系统，建立函数定义改变向量y的比率。

1. function dy=exampleode(t,y)
   
2. % function to be integrated
   
3. dy=zeros(2,1);
   
4. dy(1)=y(1)+y(2);
   
5. dy(2)=y(1); % 选择（？Alternatively）
   
6. % A=[1 1;1 0]; dy=A*y;

复制代码

利用MATLAB提供的数值求解器函数，先采用ODE45:

1. xspan = [0 10]；
   
2. ynot = [1 0];
   
3. [X,Y]=ode45(@exampleode,xspan,ynot);

复制代码

这建立了一个时间向量X(或者X的表达式)以及相应的Y向量，简单地说就是在时间X取得Y。在上面的例子中，Y的第一列为u，第二列是v。plot(X,Y):


II． 考虑这个二阶系统
第一步：通过引入向量将这个二阶ODEs变成一阶系统
第二步：将上面的那个方程系统改写成如下

将这在MATLAB中写成如下的形式：

1. function dy = secondode(x,y)
   
2. % function to be integrated
   
3. dy = zeros(4,1);
   
4. dy(1) = y(2);
   
5. dy(2) = -3*y(1)- exp(x)*y(4) + exp(2*x);
   
6. dy(3) = y(4);
   
7. dy(4) = -y(1) – cos(x)*y(2)+ sin(x);

复制代码

注意：记得将变量x改成t（它是简单独立变量）。

好，现在用ODE45解这个系统，初始状态 u(0) = 1; u'(0) = 2; v(0) = 3; v'(0) = 4；从x = 0 开始迭代到x = 3 。你需要用到下面的命令：

1. xspan = [0 3];
   
2. y0 = [1;2;3;4];
   
3. [x,y] ode45(@secondode,xpan,y0); % y的第一列是关于x的u的值，y的第二列
   
4. %是u'的值，如此类推
   
5. plot(x,y); legend('u','u’','v','v’'):

复制代码

III． 考虑下面的方程系统

其中，A,B,C与D是矩阵，y是向量。例如：

你可以降低方程的阶次，采用数值求解。首先定义

这样就可以将（1）式改写为
http://www.yculblog.com/photo/g/genial/ode9.JPG.jpg
注意：对于隐含（?implicit）求解器，例如ODE15S,ODE23T和ODE23TB，你可以将A表示成质量矩阵（？mass matrix），这是微分代数方程经常做的。
你可以写成 如下形式。
http://www.yculblog.com/photo/g/genial/ode10.JPG.jpg
你可以通过这些来链接ODE45或者另外的ODE求解器来求取数值解。例如，假设矩阵A,B,C,D如下：

可以用ODE45解这个系统如下：

1. function dy = matrixode(t,y)
   
2. % function to be integrated
   
3. dy = zeros(4,1);
   
4. dy(1) = y(3);
   
5. dy(2) = y(4);
   
6. dy(3) = -0.5*y(1)-y(2)+0.5*y(3)+y(4); % 原文此处有小错
   
7. dy(4) = -0.5*y(1)+0.5*y(3)+1;

复制代码

带有初识条件 x1(0) = 9, x2(0)=7,x1'(0)=5,x2'(0)=3,时间范围是 t=0到t=5，解决这个系统的命令如下：

1. tspan = [0 5];
   
2. x_init = [ 9; 7; 5; 3];
   
3. [t,x] = ode45(@matrixode,tspan,x_init);

复制代码

将x（2）的值用红色画出，x2'(t)的值用绿色画出，采用的命令如下：

1. plot(t,x(:,2),'r-',t,x(:,4),'g-');

复制代码

绘出图形如下页所示。

补充实例：

%------------------------------------------------------------------------
% 根据方程组写微分方程如下:
%

1. function dudt=changode(t,u)
   
2. dudt=zeros(6,1);
   
3. dudt(1)=u(4);
   
4. dudt(2)=u(5);
   
5. dudt(3)=u(6);
   
6. dudt(4)=-2*u(1)-u(2)-u(4)+3*sin(t);
   
7. dudt(5)=0.5*u(1)-u(2)+0.5*u(4)-0.5*u(6)+sin(t);
   
8. dudt(6)=(1/3)*(-4*u(3)-u(5))+2*sin(t);
   
9. dudt=dudt(;

复制代码

％――――――――――――――――――――――――――――――――――

调用ode45求解该方程组，并画出结果图：

1. function changpingpin
   
2. u0=[-1;-3;-2;1;2;3];
   
3. tspan=[0 20];
   
4. [t,u]=ode45(@changode,tspan,u0);
   
5. plot(t,u(:,1),t,u(:,2),t,u(:,3));
   
6. legend('u(1)','u(2)','u(3)');
   
7. xlabel('t');
   
8. ylabel('u');

复制代码

siyanger

这个不错，先细细看看。不过有时理论上可以，实际中不一定能行啊！

cdwxg

http://forum.vibunion.com/thread-6586-1-1.html
看下这个

zhaozbj

楼住，写得确实比较详细。可是，好像还没有写完整吧。我想看看怎么求解边值问题，特别是求解方程组的情况说明，请楼主就多传一点吧

studyboy

推荐大家使用inline函数定义微分方程（组）。呵呵。这个很不错！

songling02

用MATLAB解常微分方程的初值问题y'=-2xy,0<x<3,y(0)=1,计算解析解和数值解

hunter_009

斑主好有耐心啊，

fandalei

步长如何选择?

wangzi629

好东西，多谢分享

cheungyc

比较晦涩难懂

wei343

不错。谢谢。

octopussheng

讲的还是比较详细的，很有用处！:lol

slw

大家好！请教延时微分方程DDE23的解法，谢谢！

zhhf_3

我按照版主的指导做了一个verderpol的例子。怎么出现了问题，主要是我想让方程中的参数变化该怎么做？总是出错。
首先
function xprime = verderpol(t,x,u)
xprime = [x(2);u*(1-x(1)^2)*x(2)-x(1)];
然后在命令行中执行：
u =8;
Y0=[1;0];
options=[];

[t,x] = ode45('verderpol',[0:5:40],Y0,options,u);

如何用ode传递函数中的参数呀？

高手请帮忙。因为我要让参数不断的变化。

[ 本帖最后由 zhhf_3 于 2007-8-11 14:18 编辑 ]

vib

这些图片我怎没看不到阿？总结的很好，学习中