用MATLAB解常微分方程.docx
- 文档编号:11744284
- 上传时间:2023-03-31
- 格式:DOCX
- 页数:13
- 大小:132.55KB
用MATLAB解常微分方程.docx
《用MATLAB解常微分方程.docx》由会员分享,可在线阅读,更多相关《用MATLAB解常微分方程.docx(13页珍藏版)》请在冰豆网上搜索。
用MATLAB解常微分方程
实验四求微分方程的解
一、问题背景与实验目的
实际应用问题通过数学建模所归纳而得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法,既要研究微分方程(组)的解析解法(精确解),更要研究微分方程(组)的数值解法(近似解).
对微分方程(组)的解析解法(精确解),Matlab有专门的函数可以用,本实验将作一定的介绍.
本实验将主要研究微分方程(组)的数值解法(近似解),重点介绍Euler折线法.
二、相关函数(命令)及简介
1.dsolve('equ1','equ2';…):
Matlab求微分方程的解析解.equ1、equ2、…为方程(或条件).写方程(或条件)时用Dy表示y关于自变量的一阶导数,用用D2y表示y关于自变量的二阶导数,依此类推.
2.simplify(s):
对表达式s使用maple的化简规则进行化简.
例如:
symsx
simplify(sin(x)A2+cos(x)^2)
ans=1
3.[r,how]=simple(s):
由于Matlab提供了多种化简规则,simpie命令就是对表达式s用各种规则进行化简,然后用r返回最简形式,how返回形成这种形式所用的规则.
例如:
symsx
[r,how]=simple(cos(x)A2-sin(x)A2)
r=cos(2*x)
how=combine
4.[T,Y]=solver(odefun,tspan,y))求微分方程的数值解.
说明:
(1)其中的solver为命令ode45、ode23、ode113ode15sode23s、ode23t、ode23tb之一.
⑶在积分区间tspan=[t0,tf]上,从to到tf,用初始条件y求解.
⑷要获得问题在其他指定时间点t0,t1,t2,上的解,则令tspan=
[tott,,tf](要求是单调的).
(5)因为没有一种算法可以有效地解决所有的ODE问题,为此,Matlab提
供了多种求解器Solver,对于不同的ODE问题,采用不同的Solver.
求解器
Solver
ODE类型
特点
说明
ode45
非刚性
单步算法;4、5阶Runge-Kutta
方程;累计截断误差达(x)3
大部分场合的首选算法
ode23
非刚性
单步算法;2、3阶Runge-Kutta
方程;累计截断误差达(x)3
使用于精度较低的情形
ode113
非刚性
多步法;Adams算法;高低精
度均可到103〜106
计算时间比ode45短
ode23t
适度刚性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法;Gear's反向数值微分;精度中等
若ode45失效时,可
尝试使用
ode23s
刚性
单步法;2阶Rosebrock算法;低精度
当精度较低时,计算时
间比ode15s短
ode23tb
刚性
梯形算法;低精度
当精度较低时,计算时
间比ode15s短
⑹要特别的是:
ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的Matlab的常用程序,其中:
ode23采用龙格-库塔2阶算法,用3阶公式作误差估计来调节步长,具有低等的精度.
ode45则采用龙格-库塔4阶算法,用5阶公式作误差估计来调节步长,有中等的精度.
5.
的符
注意
ezplot(x,y,[tmin,tmax]):
符号函数的作图命令.x,y为关于参数t号函数,[tmin,tmax]为t的取值范围.
6.iniine():
建立一个内联函数.格式:
inline('expr','var1','var2',…),括号里的表达式要加引号.
例:
Q=dblquad(inline('y*sin(x)'),pi,2*pi,0,pi)
三、实验内容
2
x,并加以验证.
1.几个可以直接用Matlab求微分方程精确解的例子:
例1:
求解微分方程dy2xyxedx
求解本问题的Matlab程序为:
%line1
%line2
%line3
%line4
symsxy
y=dsolve('Dy+2*x*y=x*exp(-xA2)','x')
diff(y,x)+2*x*y-x*ex卩(帜八2)
simplify(diff(y,x)+2*x*y-x*ex卩(帜八2))
说明:
(1)行Iine1是用命令定义x,y为符号变量.这里可以不写,但为确保正确性,建议写上;
⑵行Iine2是用命令求出的微分方程的解:
1/2*exp(-xA2)*xA2+exp(-xA2)*C1
(3)行Iine3使用所求得的解.这里是将解代入原微分方程,结果应该为0,但这里给出:
⑷行Iine4用simplify()
函数对上式进行化简,结果为0,表明yy(x)的
-x^3*exp(-xA2)-2*x*exp(-xA2)*C1+2*x*(1/2*exp(-xA2)*xA2+exp(-xA2)*C1)
例2:
求微分方程xy'y
确是微分方程的解.
ex0在初始条件y
(1)2e下的特解,并画出解函
数的图形.
求解本问题的Matlab程序为:
symsxy
y=dsolve('x*Dy+y-exp(x)=0','y
(1)=2*exp
(1)','x')
ezplot(y)
eex
微分方程的特解为:
y=1/x*exp(x)+1/x*exp
(1)(Matlab格式),即卩y解函数的图形如图1:
-30
1/xexp(x)+1/xexp
(1)
50
40
30
20
10
0
-10
-20
-6-4-20246
x
2x2x的数值解,求解范围为区
0.19000.2400
0.49000.5000
并画出解函数的图形.
求解本问题的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.用ode23ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)•
例4:
求解微分方程初值问题dy22
y(0)1
间[0,0.5].
fun=iniineC-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.09000.1400
0.29000.34000.39000.4400
>>y'
ans=
1.0000
0.6440
图形结果为图
0.92470.84340.77540.71990.6764
0.6
0.95
0.9
0.85
0.8
0.75
0.7
0.65
0.15
0.20.250.30.350.40.450.5
00.05
0.62220.61050.60840.61540.6179
2.
3.用Euler折线法求解
前面讲到过,能够求解的微分方程也是十分有限的.下面介绍用Euler折线法求微分方程的数值解(近似解)的方法.
Euler折线法求解的基本思想是将微分方程初值问题
dxS
y(x0)y。
于是:
y(Xkh)y(Xk)
dy
—y
dx
y(0)1
的数值解(步长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.
求微分方程(X21)y'2xysinx0的通解.
求微分方程y''2y'5yexsinx的通解.
在初始条件xlt01,yIt00下的特解,并画出解函数yf(x)的图形.
4.
5.用Euler折线法求解微分方程初值问题
12x2
3—
y
1
y'y
分别用ode23ode45求上述第3题中的微分方程初值问题的数值解(近似解),求解区间为t[0,2].利用画图来比较两种求解器之间的差异.
y(0)
[0,2].
的数值解(步长h取0.1),求解范围为区间
6.用四阶Runge-Kutta法求解微分方程初值问题y'yexcosx,y(0)1
2.试用该方法求解第5题中的初值问题.题的常微分方程初值问题的数值解(近似解),
的数值解(步长h取0.1),求解范围为区间[0,3].
四阶
Runge-Kutta
法的迭代公式为
(Euler折线法实为一阶
yo
y(X0)
Xk1
Xk
h,
yk1
yk
6(L1
2L22L3L4)
L1
f(Xk,
yk)
k
0,1,2,,n1
L2
f(Xk
h
尹k
L3
f(xk
h
2,yk
扣)
L4
f(Xk
h,yk
hL3)
Runge-Kutta
法):
相应的Matlab程序参见附录
7.用ode45方法求上述第6从而利用画图来比较两者间的差异.
五、附录
附录1:
(fulul.m)clearf=sym('y+2*x/yA2');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))
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 微分方程