Matlab基础：实验四 求微分方程的解.doc

实验四　求微分方程的解

一、问题背景与实验目的
实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．
对微分方程（组）的解析解法(精确解)，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 题的常微分方程初值问题的数值解(近似解)，从而利用画图来比较两者间的差异．
五、附录
附录 1：(fulu1.m)
clear
f=sym('y+2*x/y^2');
a=0;
b=2;
h=0.4;
n=(b-a)/h+1;
x=0;
y=1;
szj=[x,y];
for i=1:n-1
y=y+h*subs(f,{'x','y'},{x,y});
x=x+h;
szj=[szj;x,y];
end
szj
plot(szj(:,1),szj(:,2))

附录 2：(fulu2.m)
clear
f=sym('y-exp(x)*cos(x)');
a=0;
b=3;
h=0.1;
n=(b-a)/h+1;
x=0;
y=1;
szj=[x,y];
for i=1:n-1
l1=subs(f,{'x','y'},{x,y});
l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});
l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});
l4=subs(f,{'x','y'},{x+h,y+l3*h});
y=y+h*(l1+2*l2+2*l3+l4)/6;
x=x+h;
szj=[szj;x,y];
end
szj
plot(szj(:,1),szj(:,2))
８０