实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程(组)．这就要求我们必须研究微分方程(组)的解法，既要研究微分方程(组)的解析解法(精确解)，更要研究微分方程(组)的数值解法(近似解)．

对微分方程(组)的解析解法(精确解)，Matlab有专门的函数可以用，本实验将作一定的介绍．

本实验将主要研究微分方程(组)的数值解法(近似解)，重点介绍Euler折线法．

1．dsolve('equ1','equ2',…)：Matlab求微分方程的解析解．equ1、equ2、…为方程(或条件)．写方程(或条件)时用Dy表示y关于自变量的一阶导数，用用D2y表示y关于自变量的二阶导数，依此类推．

2．simplify(s)：对表达式s使用maple的化简规则进行化简．

例如：

syms x

simplify(sin(x)^2 + cos(x)^2)

ans=1

3．[r,how]=simple(s)：由于Matlab提供了多种化简规则，simple命令就是对表达式s用各种规则进行化简，然后用r返回最简形式，how返回形成这种形式所用的规则．

例如：

syms x

[r,how]=simple(cos(x)^2-sin(x)^2)

r = cos(2*x)

how = combine

4．[T,Y] = solver(odefun,tspan,y0)求微分方程的数值解．

说明：

(1)其中的solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb之一．

(2) odefun是显式常微分方程：

(3)在积分区间tspan=

上，从到，用初始条件求解．

(4)要获得问题在其他指定时间点上的解，则令tspan=

(要求是单调的)．

(5)因为没有一种算法可以有效地解决所有的ODE问题，为此，Matlab提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver．

求解器

Solver

ODE类型

特点

说明

ode45

非刚性

单步算法；4、5阶Runge-Kutta方程；累计截断误差达

大部分场合的首选算法

ode23

非刚性

单步算法；2、3阶Runge-Kutta方程；累计截断误差达

使用于精度较低的情形

ode113

非刚性

多步法；Adams算法；高低精度均可到

计算时间比ode45短

ode23t

适度刚性

采用梯形算法

适度刚性情形

ode15s

刚性

多步法；Gear's反向数值微分；精度中等

若ode45失效时，可尝试使用

ode23s

刚性

单步法；2阶Rosebrock算法；低精度

当精度较低时，计算时间比ode15s短

ode23tb

刚性

梯形算法；低精度

当精度较低时，计算时间比ode15s短

(6)要特别的是：ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的Matlab的常用程序，其中：

ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度．

ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度．

5．ezplot(x,y,[tmin,tmax])：符号函数的作图命令．x,y为关于参数t的符号函数，[tmin,tmax]为t的取值范围．

6．inline()：建立一个内联函数．格式：inline('expr', 'var1', 'var2',…)，注意括号里的表达式要加引号．

例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi).

1.几个可以直接用Matlab求微分方程精确解的例子：

例1：求解微分方程，并加以验证．

求解本问题的Matlab程序为：

syms x y                                  %line1

y=dsolve('Dy+2*x*y=x*exp(-x^2)','x')          %line2

diff(y,x)+2*x*y-x*exp(-x^2)                  %line3

simplify(diff(y,x)+2*x*y-x*exp(-x^2))          %line4

说明：

(1)行line1是用命令定义x,y为符号变量．这里可以不写，但为确保正确性，建议写上；

(2)行line2是用命令求出的微分方程的解：

1/2*exp(-x^2)*x^2+exp(-x^2)*C1

(3)行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：

-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)

(4)行line4用simplify()函数对上式进行化简，结果为0，表明的确是微分方程的解．

例2：求微分方程在初始条件下的特解，并画出解函数的图形．

求解本问题的Matlab程序为：

syms x y

y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')

ezplot(y)

微分方程的特解为：y=1/x*exp(x)+1/x* exp (1)(Matlab格式)，即，解函数的图形如图1：

图1

例3：求微分方程组在初始条件下的特解，并画出解函数的图形．

求解本问题的Matlab程序为：

syms x y t

[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')

simple(x);

simple(y);

ezplot(x,y,[0,1.3]);axis auto

微分方程的特解(式子特别长)以及解函数的图形均略．

2.用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．

例4：求解微分方程初值问题的数值解，求解范围为区间[0, 0.5]．

fun=inline('-2*y+2*x^2+2*x','x','y');

[x,y]=ode23(fun,[0,0.5],1);

x';

y';

plot(x,y,'o-')

>> x'

ans =

0.0000   0.0400   0.0900   0.1400   0.1900   0.2400

0.2900   0.3400   0.3900   0.4400   0.4900   0.5000

>> y'

ans =

1.0000   0.9247   0.8434   0.7754   0.7199   0.6764

0.6440   0.6222   0.6105   0.6084   0.6154   0.6179

图形结果为图2．

图2

例5：求解描述振荡器的经典的Ver der Pol微分方程

分析：令则

先编写函数文件verderpol.m：

function xprime = verderpol(t,x)

global mu;

xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)];

再编写命令文件vdp1.m：

global mu;

mu = 7;

y0=[1;0]

[t,x] = ode45('verderpol',[0,40],y0);

x1=x(:,1);x2=x(:,2);

plot(t,x1)

图形结果为图3．

图3

3.用Euler折线法求解

前面讲到过，能够求解的微分方程也是十分有限的．下面介绍用Euler折线法求微分方程的数值解(近似解)的方法．

Euler折线法求解的基本思想是将微分方程初值问题

化成一个代数方程，即差分方程，主要步骤是用差商替代微商,于是：

记，从而，则有

例6：用Euler折线法求解微分方程初值问题

的数值解(步长h取0.4)，求解范围为区间[0,2]．

解：本问题的差分方程为

相应的Matlab程序见附录1．

数据结果为：

0         1.0000

0.4000    1.4000

0.8000    2.1233

1.2000    3.1145

1.6000    4.4593

2.0000    6.3074

图形结果见图4：

图4

特别说明：本问题可进一步利用四阶Runge-Kutta法求解，读者可将两个结果在一个图中显示，并和精确值比较，看看哪个更“精确”？(相应的Matlab程序参见附录2)．

1.求微分方程的通解．

2.求微分方程的通解．

3.求微分方程组

在初始条件下的特解，并画出解函数的图形．

4.分别用ode23、ode45求上述第3题中的微分方程初值问题的数值解(近似解)，求解区间为．利用画图来比较两种求解器之间的差异．

5.用Euler折线法求解微分方程初值问题

的数值解(步长h取0.1)，求解范围为区间[0,2]．

6.用四阶Runge-Kutta法求解微分方程初值问题

的数值解(步长h取0.1)，求解范围为区间[0,3]．

四阶Runge-Kutta法的迭代公式为(Euler折线法实为一阶Runge-Kutta法)：

相应的Matlab程序参见附录2．试用该方法求解第5题中的初值问题．

7.用ode45方法求上述第6题的常微分方程初值问题的数值解(近似解)，从而利用画图来比较两者间的差异．