实验四 求微分方程的解
一、问题背景与实验目的
实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．
对微分方程（组）的解析解法(精确解)，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,y 0) 求微分方程的数值解． 说明：
(1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一．
(2) odefun 是显式常微分方程：⎪⎩⎪⎨⎧==0
0)()
,(y t y y t f dt dy
(3) 在积分区间 tspan =],[0f t t 上，从0t 到f t ，用初始条件0y 求解．
(4) 要获得问题在其他指定时间点 ,210,,t t t 上的解，则令 tspan =
],,,[,210f t t t t （要求是单调的）．
(5) 因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver ，对于不同的ODE 问题，采用不同的Solver ．
(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：求解微分方程22x xe xy dx
dy
-=+，并加以验证．
求解本问题的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， 表明
)(x y y =的确是微分方程的解．
例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．
求解本问题的 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 格式)，即
x
e e y x +=，解函数的图形如图 1：
图1
例3：求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++035y x dt
dy e y x dt
dx t
在初始条件0|,1|00====t t y x 下的特解，
并画出解函数的图形．
求解本问题的 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：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222
y x
x y dx dy 的数值解，求解范围为区
间[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
