用MATLAB解常微分方程.docx
- 文档编号:26216589
- 上传时间:2023-06-17
- 格式:DOCX
- 页数:12
- 大小:161.30KB
用MATLAB解常微分方程.docx
《用MATLAB解常微分方程.docx》由会员分享,可在线阅读,更多相关《用MATLAB解常微分方程.docx(12页珍藏版)》请在冰豆网上搜索。
用MATLAB解常微分方程
用MATLAB解常微分方程
LT
三、实验内容
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.09000.14000.19000.2400
0.29000.34000.39000.44000.49000.5000
>>y'
ans=
1.00000.92470.84340.77540.71990.6764
0.64400.62220.61050.60840.61540.6179
图形结果为图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
0.40001.4000
0.80002.1233
1.20003.1145
1.60004.4593
2.00006.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];
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 微分方程