(精品)第6章微分方程问题的解法.ppt
《(精品)第6章微分方程问题的解法.ppt》由会员分享,可在线阅读,更多相关《(精品)第6章微分方程问题的解法.ppt(86页珍藏版)》请在得力文库 - 分享文档赚钱的网站上搜索。
1、第六章 微分方程问题的解法微分方程的解析解方法常微分方程问题的数值解法微分方程问题算法概述四阶定步长 Runge-Kutta算法及 MATLAB 实现一阶微分方程组的数值解微分方程转换特殊微分方程的数值解边值问题的计算机求解偏微分方程的解6.1 微分方程的解析解方法格式:y=dsolve(f1,f2,fm)格式:指明自变量 y=dsolve(f1,f2,fm,x)fi即可以描述微分方程,又可描述初始条件或边界条件。如:描述微分方程时 描述条件时 例:syms t;u=exp(-5*t)*cos(2*t+1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*uuu=87*exp
2、(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10 syms t y;y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=,.87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10)y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=,.87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1).+10,y(0)=3,Dy(0)=2,D2y(0)=0,D3y(0)=0)分别处理系数,如:n,d=rat(double(vpa(-4
3、45/26*cos(1)-51/13*sin(1)-69/2)ans=-8704 185%rat()最接近有理数的分数判断误差:vpa(-445/26*cos(sym(1)-51/13*sin(1)-69/2+8704/185)ans=.114731975864790922564144636e-4 y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=,.87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+.10,y(0)=1/2,Dy(pi)=1,D2y(2*pi)=0,Dy(2*pi)=1/5);如果用推导的方法求Ci的值,
4、每个系数的解析解至少要写出10数行,故可采用有理式近似 的方式表示.vpa(y,10)%有理近似值ans=1.196361839*exp(-5.*t)+.4166666667-.4785447354*sin(t)*cos(t)*exp(-5.*t)-.4519262218e-1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)2*exp(-5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(-3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(-1.*t)+442590
5、.9059*exp(-4.*t)例:求解 x,y=dsolve(D2x+2*Dx=x+2*y-exp(-t),Dy=4*x+3*y+4*exp(-t)例:syms t x x=dsolve(Dx=x*(1-x2)x=1/(1+exp(-2*t)*C1)(1/2)-1/(1+exp(-2*t)*C1)(1/2)syms t x;x=dsolve(Dx=x*(1-x2)+1)Warning:Explicit solution could not be found;implicit solution returned.In D:MATLAB6p5toolboxsymbolicdsolve.m at
6、line 292x=t-Int(1/(a-a3+1),a=.x)+C1=0故只有部分非线性微分方程有解析解。6.2 微分方程问题的数值解法6.2.1 微分方程问题算法概述微分方程求解的误差与步长问题:6.2.2 四阶定步长Runge-Kutta算法 及 MATLAB 实现 function tout,yout=rk_4(odefile,tspan,y0)y0初值列向量 t0=tspan(1);th=tspan(2);if length(tspan)t_final=100;x0=0;0;1e-10;%t_final为设定的仿真终止时间 t,x=ode45(lorenzeq,0,t_final,x
7、0);plot(t,x),figure;%打开新图形窗口 plot3(x(:,1),x(:,2),x(:,3);axis(10 42-20 20-20 25);%根据实际数值手动设置坐标系可采用comet3()函数绘制动画式的轨迹。comet3(x(:,1),x(:,2),x(:,3)描述微分方程是常微分方程初值问题数值求解的关键。f1=inline(-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);,.-x(1)*x(2)+28*x(2)-x(3),t,x);t_final=100;x0=0;0;1e-10;t,x=ode45(f1,0,t_final,x0);plo
8、t(t,x),figure;plot3(x(:,1),x(:,2),x(:,3);axis(10 42-20 20-20 25);得出完全一致的结果。6.2.3.3 MATLAB 下带有附加参数的微分方程求解例:编写函数function xdot=lorenz1(t,x,flag,beta,rho,sigma)flag变量是不能省略的 xdot=-beta*x(1)+x(2)*x(3);-rho*x(2)+rho*x(3);-x(1)*x(2)+sigma*x(2)-x(3);求微分方程:t_final=100;x0=0;0;1e-10;b2=2;r2=5;s2=20;t2,x2=ode45(
9、lorenz1,0,t_final,x0,b2,r2,s2);plot(t2,x2),options位置为,表示不需修改控制选项 figure;plot3(x2(:,1),x2(:,2),x2(:,3);axis(0 72-20 22 -35 40);f2=inline(-beta*x(1)+x(2)*x(3);-rho*x(2)+rho*x(3);,.-x(1)*x(2)+sigma*x(2)-x(3),t,x,flag,beta,rho,sigma);flag变量是不能省略的6.2.4 微分方程转换6.2.4.1 单个高阶常微分方程处理方法例:函数描述为:function y=vdp_eq
10、(t,x,flag,mu)y=x(2);-mu*(x(1)2-1)*x(2)-x(1);x0=-0.2;-0.7;t_final=20;mu=1;t1,y1=ode45(vdp_eq,0,t_final,x0,mu);mu=2;t2,y2=ode45(vdp_eq,0,t_final,x0,mu);plot(t1,y1,t2,y2,:)figure;plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:)x0=2;0;t_final=3000;mu=1000;t,y=ode45(vdp_eq,0,t_final,x0,mu);由于变步长所采用的步长过小,所需时间较长,导致
11、输出的y矩阵过大,超出计算机存储空间容量。所以不适合采用ode45()来求解,可用刚性方程求解算法ode15s()。6.2.4.2 高阶常微分方程组的变换方法例:描述函数:function dx=apolloeq(t,x)mu=1/82.45;mu1=1-mu;r1=sqrt(x(1)+mu)2+x(3)2);r2=sqrt(x(1)-mu1)2+x(3)2);dx=x(2);2*x(4)+x(1)-mu1*(x(1)+mu)/r13-mu*(x(1)-mu1)/r23;x(4);-2*x(2)+x(3)-mu1*x(3)/r13-mu*x(3)/r23;求解:x0=1.2;0;0;-1.04
12、935751;tic,t,y=ode45(apolloeq,0,20,x0);tocelapsed_time=0.8310 length(t),plot(y(:,1),y(:,3)ans=689得出的轨道不正确,默认精度RelTol设置得太大,从而导致的误差传递,可减小该值。改变精度:options=odeset;options.RelTol=1e-6;tic,t1,y1=ode45(apolloeq,0,20,x0,options);tocelapsed_time=0.8110 length(t1),plot(y1(:,1),y1(:,3),ans=1873 min(diff(t1)ans=
13、1.8927e-004 plot(t1(1:end-1),diff(t1)例:x0=1.2;0;0;-1.04935751;tic,t1,y1=rk_4(apolloeq,0,20,0.01,x0);tocelapsed_time=4.2570 plot(y1(:,1),y1(:,3)%绘制出轨迹曲线显而易见,这样求解是错误的,应该采用更小的步长。tic,t2,y2=rk_4(apolloeq,0,20,0.001,x0);tocelapsed_time=124.4990 计算时间过长 plot(y2(:,1),y2(:,3)%绘制出轨迹曲线严格说来某些点仍不满足106的误差限,所以求解常微分
14、方程组时建议采用变步长算法,而不是定步长算法。例:用MATLAB符号工具箱求解,令 syms x1 x2 x3 x4 dx,dy=solve(dx+2*x4*x1=2*dy,dx*x4+3*x2*dy+x1*x4-x3=5,dx,dy)%dx,dy为指定变量dx=-2*(3*x4*x1*x2+x4*x1-x3-5)/(2*x4+3*x2)dy=(2*x42*x1-x4*x1+x3+5)/(2*x4+3*x2)对于更复杂的问题来说,手工变换的难度将很大,所以如有可能,可采用计算机去求解有关方程,获得解析解。如不能得到解析解,也需要在描写一阶常微分方程组时列写出式子,得出问题的数值解。6.3特殊微
15、分方程的数值解 6.3.1 刚性微分方程的求解刚性微分方程 一类特殊的常微分方程,其中一些解变化缓慢,另一些变化快,且相差悬殊,这类方程常常称为刚性方程。MATLAB采用求解函数ode15s(),该函数的调用格式和ode45()完全一致。t,x=ode15s(Fun,t0,tf,x0,options,p1,p2,)例:计算 h_opt=odeset;h_opt.RelTol=1e-6;x0=2;0;t_final=3000;tic,mu=1000;t,y=ode15s(vdp_eq,0,t_final,x0,h_opt,mu);tocelapsed_time=2.5240作图 plot(t,y
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 精品 微分方程 问题 解法
限制150内