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*d
2、iff(u,t)+2*u uu = 87*exp(-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)
3、=0),分别处理系数，如： n,d=rat(double(vpa(-445/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
4、(pi)=1,D2y(2*pi)=0,Dy(2*pi)=1/5);如果用推导的方法求Ci的值，每个系数的解析解至少要写出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)+
5、31319.63786*exp(-2.*t)-219.1293619*exp(-1.*t)+442590.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; imp
6、licit solution returned. In D:MATLAB6p5toolboxsymbolicdsolve.m at line 292 x = 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)
7、=3, h=tspan(3); % tspan=t0,th,helse, h=tspan(2)-tspan(1); th=tspan(end); end 等间距数组 tout=t0:h:th; yout=;for t=toutk1=h*eval(odefile (t,y0); % odefile是一个字符串变量，为表示微分方程f( )的文件名。k2=h*eval(odefile (t+h/2,y0+0.5*k1);k3=h*eval(odefile (t+h/2,y0+0.5*k2);k4=h*eval(odefile (t+h,y0+k3);y0=y0+(k1+2*k2+2*k3+k4)/6
8、;yout=yout; y0;end 由效果看，该算法不是一个较好的方法。,6.2.3 一阶微分方程组的数值解 6.2.3.1 四阶五级Runge-Kutta-Felhberg算法,通过误差向量调节步长，此为自动变步长方法。 四阶五级RKF算法有参量系数表。,6.2.3.2 基于 MATLAB 的微分方程,求解函数 格式1： 直接求解 t,x=ode45(Fun,t0,tf,x0) 格式2： 带有控制参数 t,x=ode45(Fun,t0,tf,x0，options) 格式3： 带有附加参数 t,x=ode45(Fun,t0,tf,x0,options,p1,p2,) t0,tf求解区间， x
9、0初值问题的初始状态变量。,描述需要求解的微分方程组：不需附加变量的格式function xd=funname(t,x)可以使用附加变量function xd=funname(t,x,flag,p1,p2,) t是时间变量或自变量（必须给），x为状态向量， xd为返回状态向量的导数。flag用来控制求解过程,指定初值，即使初值不用指定，也必须有该变量占位。 修改变量：options唯一结构体变量，用odeset( )修改。options=odeset(RelTol,1e-7);options= odeset; options. RelTol= 1e-7;,例：自变函数function xdot
10、 = lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3);, t_final=100; x0=0;0;1e-10; % t_final为设定的仿真终止时间 t,x=ode45(lorenzeq,0,t_final,x0); plot(t,x), figure; % 打开新图形窗口 plot3(x(:,1),x(:,2),x(:,3); axis(10 42 -20 20 -20 25); % 根据实际数值手动设置坐标系,可采用comet3( )函数绘制动画式的轨迹。 comet3(x(:
11、,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); plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3); axis(10 42 -20 20 -20 25);得出完全一致的结果。,6.2.3.3 MATLAB 下带有附加参数的微分方程求解,例：,编写函数
12、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(lorenz1,0,t_final,x0,b2,r2,s2); plot(t2,x2), options位置为，表示不需修改控制选项 figure; plot3(x2(:,1),x2(:,2),x2(:,3
13、); 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(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(
14、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);由于变步长所采用的步长过小，所需时间较长，导致输出的y矩阵过大，超出计算机存储空间容量。所以不适合采用ode45()来求解，可用刚性方程求解算法ode15s( )。,6.2.4.2 高阶常微分方程
15、组的变换方法,例：,描述函数：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.04935751; tic, t,y=ode45(apolloeq,0,20,x0); toc elapsed_time =0.83
16、10 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); toc elapsed_time =0.8110 length(t1), plot(y1(:,1),y1(:,3), ans =1873, min(diff(t1) ans =1.8927e-004 plot(t1(1:end-1),diff(t1)
17、,例：, x0=1.2; 0; 0; -1.04935751; tic, t1,y1=rk_4(apolloeq,0,20,0.01,x0); toc elapsed_time =4.2570 plot(y1(:,1),y1(:,3) % 绘制出轨迹曲线显而易见，这样求解 是错误的，应该采用 更小的步长。, tic, t2,y2=rk_4(apolloeq,0,20,0.001,x0); toc elapsed_time =124.4990 计算时间过长 plot(y2(:,1),y2(:,3) % 绘制出轨迹曲线严格说来某些点仍不 满足106的误差限， 所以求解常微分方程 组时建议采用变步长
18、 算法，而不是定步长 算法。,例：,用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
19、特殊微分方程的数值解 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); toc elapsed_time =2.52
20、40,作图 plot(t,y(:,1); figure; plot(t,y(:,2)y(:,1)曲线变化较平滑, y(:,2)变化在某些点上较快。,例：定义函数 function dy=c7exstf2(t,y)dy=0.04*(1-y(1)-(1-y(2)*y(1)+0.0001*(1-y(2)2;-104*y(1)+3000*(1-y(2)2;,方法一 tic,t2,y2=ode45(c7exstf2,0,100,0;1); toc elapsed_time =229.4700 length(t2), plot(t2,y2) ans =356941,步长分析： format long, m
21、in(diff(t2), max(diff(t2) ans =0.00022220693884 0.00214971787184 plot(t2(1:end-1),diff(t2),方法二，用ode15s()代替ode45() opt=odeset; opt.RelTol=1e-6; tic,t1,y1=ode15s(c7exstf2,0,100,0;1,opt); toc elapsed_time =0.49100000000000 length(t1), plot(t1,y1) ans =169,6.3.2 隐式微分方程求解,隐式微分方程为不能转化为显式常微分方程组的方程 例：,编写函数：
22、 function dx=c7ximp(t,x)A=sin(x(1) cos(x(2); -cos(x(2) sin(x(1);B=1-x(1); -x(2); dx=inv(A)*B; 求解： opt=odeset; opt.RelTol=1e-6; t,x=ode45(c7ximp,0,10,0; 0,opt); plot(t,x),6.3.3 微分代数方程求解,例：,编写函数function dx=c7eqdae(t,x)dx=-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);x(1)+x(2)+x
23、(3)-1; M=1,0,0; 0,1,0; 0,0,0; options=odeset; options.Mass=M; Mass微分代数方程中 的质量矩阵（控制参数） x0=0.8; 0.1; 0.1; t,x=ode15s(c7eqdae,0,20,x0,options); plot(t,x),编写函数： function dx=c7eqdae1(t,x)dx=-0.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2);2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2);, x0=0.8; 0.1; fDae=inline(-0
24、.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2);,. 2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2),t,x); t1,x1=ode45(fDae,0,20,x0); plot(t1,x1,t1,1-sum(x1),6.3.3延迟微分方程求解,sol:结构体数据，sol.x:时间向量t, sol.y:状态向量。,例：,编写函数： function dx=c7exdde(t,x,z)xlag1=z(:,1); %第一列表示提取xlag2=z(:,2);dx=1-3*x(1)-xlag1(2)-0.2*xlag2(1)3-xl
25、ag2(1);x(3); 4*x(1)-2*x(2)-3*x(3);历史数据函数： function S=c7exhist(t)S=zeros(3,1);,求解： lags=1 0.5; tx=dde23(c7exdde,lags,zeros(3,1),0,10); plot(tx.x,tx.y(2,:) 与ode45()等返回的x矩阵不一样，它是按行排列的。,6.4边值问题的计算机求解,6.4.1 边值问题的打靶算法,数学方法描述：以二阶方程为例,编写函数： 线性的 function t,y=shooting(f1,f2,tspan,x0f,varargin)t0=tspan(1); tfi
26、nal=tspan(2); ga=x0f(1); gb=x0f(2);t,y1=ode45(f1,tspan,1;0,varargin); t,y2=ode45(f1,tspan,0;1,varargin);t,yp=ode45(f2,tspan,0;0,varargin); m=(gb-ga*y1(end,1)-yp(end,1)/y2(end,1);t,y=ode45(f2,tspan,ga;m,varargin);,例：编写函数： function xdot=c7fun1(t,x)xdot=x(2); -2*x(1)+3*x(2);function xdot=c7fun2(t,x)xdo
27、t=x(2); t-2*x(1)+3*x(2); t,y=shooting(c7fun1, c7fun2,0,1,1;2); plot(t,y),原方程的解析解为解的检验 y0=(exp(2)-3)*exp(t)+(3-exp(1)*exp(2*t)/(4*exp(1)*(exp(1)-1)+3/4+t/2; norm(y(:,1)-y0) % 整个解函数检验 ans =4.4790e-008 norm(y(end,1)-2) % 终点条件检验 ans =2.2620e-008,非线性方程边值问题的打靶算法：用Newton迭代法处理,编写函数： function t,y=nlbound(fun
28、cn,funcv,tspan,x0f,tol,varargin)t0=tspan(1);tfinal=tspan(2); ga=x0f(1); gb=x0f(2); m=1; m0=0;while (norm(m-m0)tol), m0=m;t,v=ode45(funcv,tspan,ga;m;0;1,varargin);m=m0-(v(end,1)-gb)/(v(end,3);endt,y=ode45(funcn,tspan,ga;m,varargin);,例：编写两个函数： function xdot=c7fun3(t,x)xdot=x(2); 2*x(1)*x(2); x(4); 2*x
29、(2)*x(3)+2*x(1)*x(4);function xdot=c7fun4(t,x)xdot=x(2); 2*x(1)*x(2);, t,y=nlbound(c7fun4,c7fun3,0,pi/2,-1,1,1e-8); plot(t,y); set(gca,xlim,0,pi/2); 精确解：检验： y0=tan(t-pi/4); norm(y(:,1)-y0) ans =1.6629e-005 norm(y(end,1)-1) ans =5.2815e-006,6.4.2 线性微分方程的有限差分算法,把等式左边用差商表示。,编写函数： function x,y=fdiff(fun
30、cs,tspan,x0f,n)t0=tspan(1);tfinal=tspan(2); ga=x0f(1); gb=x0f(2);h=(tfinal-t0)/n;for i=1:n, x(i)=t0+h*(i-1); end,x0=x(1:n-1); t=-2+h2*feval(funcs,x0,2); tmp=feval(funcs,x0,1);v=1+h*tmp/2; w=1-h*tmp/2; b=h2*feval(funcs,x0,3);b(1)=b(1)-w(1)*ga; b(n-1)=b(n-1)-v(n-1)*gb; b=b; A=diag(t);for i=1:n-2, A(i,
31、i+1)=v(i); A(i+1,i)=w(i+1); endy=inv(A)*b; x=x tfinal; y=ga; y; gb;,例：编写函数： function y=c7fun5(x,key)switch keycase 1, y=1+x;case 2, y=1-x;otherwise, y=1+x.2;end t,y=fdiff(c7fun5,0,1,1,4,50); plot(t,y),6.5 偏微分方程求解入门 6.5.1 偏微分方程组求解,函数描述：,边界条件的函数描述：初值条件的函数描述：u0=pdeic(x),例：,函数描述：,function c,f,s=c7mpde(x
32、,t,u,du)c=1;1; y=u(1)-u(2); F=exp(5.73*y)-exp(-11.46*y); s=F*-1; 1;f=0.024*du(1); 0.17*du(2);,描述边界条件的函数 function pa,qa,pb,qb=c7mpbc(xa,ua,xb,ub,t)pa=0; ua(2); qa=1;0; pb=ub(1)-1; 0; qb=0;1;,描述初值： function u0=c7mpic(x)u0=1; 0; 求解： x=0:.05:1; t=0:0.05:2; m=0; sol=pdepe(m,c7mpde,c7mpic,c7mpbc,x,t); sur
33、f(x,t,sol(:,:,1)， figure; surf(x,t,sol(:,:,2),6.5.2 二阶偏微分方程的数学描述,椭圆型偏微分方程：,抛物线型偏微分方程：双曲型偏微分方程：特征值型偏微分方程：,6.5.3 偏微分方程的求解界面应用简介 6.5.3.1 偏微分方程求解程序概述,启动偏微分方程求解界面 在 MATLAB 下键入 pdetool 该界面分为四个部分 菜单系统 工具栏 集合编辑 求解区域,6.5.3.2 偏微分方程求解区域绘制,1）用工具栏中的椭圆、矩形等绘制一些区域。 2）在集合编辑栏中修改其内容。如（R1E1E2）E3 3）单击工具栏中 按纽可得求解边界。 4）选择
34、Boundary-Remove All Subdomain Borders菜单项，消除相邻区域中间的分隔线。 5）单击 按纽可将求解区域用三角形划分成网格。可用 按纽加密。,6.5.3.3 偏微分方程边界条件描述,选择Boundary-Specify Boundary Conditions菜单,6.5.3.4 偏微分方程求解举例,例：求解： 1）绘制求解区域。 2）描述边界条件（Boundary-Specify Boundary Conditions）。 3）选择偏微分方程的类型。单击工具栏中的PDE图标，在打开的新窗口选择Hyperbolic选项，输入参数c,a,f,d. 4）求解。单击工具
35、栏中的等号按钮。,显示： 1）图形颜色表示t=0时u(x,y)的函数值。 2）单击工具栏中的三维图标将打开一新的对话框，若再选择Contour可绘制等值线，若选择Arrows选项可绘制引力线。若单独选择Height（3d-plot)，则在另一窗口绘制出三维图形。 3）可在单击三维图标打开的新对话框中，对Property栏目的各个项目重新选择。 4）可修改微分方程的边界条件，重新求解。 动画： 1）Solve-Parameters对话框时间向量改为0:0.1:2。 2）三维图标打开的对话框中选择Animation选项，单击Options按纽设置播放速度。Plot-Export Movie 菜单。,6.5.3.5 函数参数的偏微分方程求解,例：(椭圆型）,求解： 1）求解区域不变。 2）描述边界条件，u=0。 3）选择偏微分方程的类型。单击工具栏中的PDE图标，在打开的新窗口选择Elliptic选项，输入参数c=1./sqrt(1+ux.2+uy.2), a=x.2+y.2 , f=exp(-x.2-y.2). 4)再打开Solve-Parameters对话框，选定Use nonlinear solve属性（该属性只适于椭圆性偏微分方程） 5）求解。单击工具栏中的等号按钮。,
