exp求微分方程的解.docx
《exp求微分方程的解.docx》由会员分享，可在线阅读，更多相关《exp求微分方程的解.docx（14页珍藏版）》请在冰豆网上搜索。
exp求微分方程的解
exp4求微分方程的解
————————————————————————————————作者：
————————————————————————————————日期：
实验四 求微分方程的解
一、问题背景与实验目的
实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法,既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．
对微分方程（组）的解析解法（精确解），Matlab有专门的函数可以用，本实验将作一定的介绍．
本实验将主要研究微分方程（组）的数值解法（近似解）,重点介绍Euler折线法．
二、相关函数（命令）及简介
1．dsolve（’equ1'，’equ2',…）:
Matlab求微分方程的解析解．equ1、equ2、…为方程（或条件）．写方程（或条件）时用Dy表示y关于自变量的一阶导数，用用D2y表示y关于自变量的二阶导数，依此类推．
2．simplify（s）：
对表达式s使用maple的化简规则进行化简．
例如:
symsx
simplify（sin（x）^2+cos（x）^2）
ans=1
3．[r，how]=simple（s）：
由于Matlab提供了多种化简规则，simple命令就是对表达式s用各种规则进行化简，然后用r返回最简形式,how返回形成这种形式所用的规则．
例如：
symsx
［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程序为：
symsxy%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程序为：
symsxy
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程序为：
symsxyt
[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]）;axisauto
微分方程的特解（式子特别长）以及解函数的图形均略．
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.00000.04000。
14000.1900
0。
29000.34000。
39000。
44000。
49000。
5000
〉〉y’
ans=
1。
00000。
92470。
84340。
77540.7199
0。
64400。
62220.61050。
60840。
6154
图形结果为图2．
图2
例5：
求解描述振荡器的经典的VerderPol微分方程
分析：
令
则
先编写函数文件verderpol。
m：
functionxprime=verderpol（t，x）
globalmu;
xprime=[x
（2）；mu＊（1-x
（1）^2）＊x
（2）-x
（1）]；
再编写命令文件vdp1。
m:
globalmu；
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．
数据结果为：
01。
0000
1.60004。
4593
图形结果见图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］；
fori=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］；
fori=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））
升级会员 