1、第五章 微分方程问题的解法,常系数线性微分方程的解析解方法 常微分方程问题的数值解法 微分方程问题算法概述 四阶定步长 Runge-Kutta算法及 MATLAB 实现 一阶微分方程组的数值解 微分方程转换 特殊微分方程的数值解 边值问题的计算机求解 偏微分方程的解,5.1 常系数线性微分方程的解析解方法,数学描述:,特征方程(多项式代数方程):,该代数方程的根称为原常系数方程的特征根 根据特征根的情况可求得原方程的解析解,格式:y=dsolve(f1, f2, , fm)格式:指明自变量y=dsolve(f1, f2, , fm ,x)fi 即可以描述微分方程,又可描述初始条件或边界条件。如
2、:描述微分方程时描述条件时,例:, syms t; u=exp(-5*t)*cos(2*t+1)+5; uu=5*diff(u,t,2)+4*diff(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)*co
3、s(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(-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+2
4、4*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的值,每个系数的解析解至少要写出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(
5、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.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
6、=dsolve(Dx=x*(1-x2)+1) Warning: Explicit solution could not be found; implicit solution returned. In D:MATLAB6p5toolboxsymbolicdsolve.m at line 292 x = t-Int(1/(a-a3+1),a=x)+C1=0 故只有部分非线性微分方程可直接求出其解析解。,5.2 微分方程问题的数值解法,5.2.1 算法概述,微分方程求解的误差与步长问题:,function outx,outy=MyEuler(fun,x0,xt,y0,PointNum)% fun
7、表示f(x,y); x0,xt:自变量的初值和终值; %y0:函数在x0处的值,其可以为向量形式; %PointNum表示自变量在x0,xt上取的点数if nargin=4 | PointNum=0 %PointNum 默认值为100PointNum=100; end if nargin=3 % y0默认值为0 y0=0; PointNum=100; end,h=(xt-x0)/PointNum; %计算步长h x=x0+0:PointNum*h; %自变量数组 y= zeros(size(x);与自变量点相应的函数值 y(1)=y0; %初值 for k = 1:PointNumy(k+1)
8、= y(k)+h*feval(fun,x(k),y(k);%对于所取的点x迭代计算y值 end outy=y; outx=x; plot(x,y) %画出方程解的函数图,例:求,f1=(x,y)sin(x)+y; 2*pi/0.4 % 计算 xt=2pi 时应取得点数 ans =15.7080 x1,y1=MyEuler(f1,0,2*pi,1,16); %h=0.4 欧拉法所的的解 x11,y11=MyEuler(f1,0,2*pi,1,32); % h=0.2 欧拉法所的的解,h=0.2和h=0.4的数值解。,y=dsolve(Dy=y+sin(t),y(0)=1); %该常微分方程的解析
9、解 for k=1:33 t(k)=x11(k); y2(k)=subs(y,t(k); %求其对应点的离散解 end plot(x1,y1,+b,x11,y11,og,x11,y2,*r) legend(h=0.4的欧拉法解,h=0.2的欧拉解,符号解),误差判断:,改进的Euler算法,改进Euler迭代公式,function Xout,Yout=MyEulerPro(fun,x0,xt,y0,PointNumber)%MyEulerPro 用改进的欧拉法解微分方程 if nargin=4 | PointNumber=0 %PointNumer默认值为100PointNumer=100;
10、end if nargin=3 %y0默认值为0y0=0; PointNumer=100; end,算法实现:,h=(xt-x0)/PointNumber; %计算所取的两离 散点之间的距离 x=x0+0:PointNumber*h; %表示出离散的自变量x,初始点为x0 y=zeros(size(x); 输入与x相应的函数值向量 y(1)=y0;初始值 for k=1:PointNumber %迭代计算过程z1=y(k)+h*feval(fun,x(k),y(k); z2=y(k)+h*feval(fun,x(k+1),z1); y(k+1)=1/2*(z1+z2); end Xout=x;
11、Yout=y;,例:对方程 应用改进欧拉法、简单欧拉方法 以及微分方程符号求解方法分别求解并比较结果; f1=(x,y)sin(x)+y; x3,y3=MyEulerPro(f1,0,2*pi,1,128); x,y1=MyEuler(f1,0,2*pi,1,128);%欧拉法所的的解 y=dsolve(Dy=y+sin(t),y(0)=1);%该常微分方程的符号解 for k=1:129 t(k)=x(k); y2(k)=subs(y,t(k);%求解析解对应点的离散解 end plot(x,y1,-b,x3,y3,og,x,y2,*r) legend(简单欧拉法解,改进欧拉解,符号解),5
12、.2.2 四阶定步长Runge-Kutta算法及Matlab实现,Euler算法:一阶Taylor多项式近似函数,基本思想:,RK一般算法:,推广:高级Taylor多项式近似函数提高数值解精度,例如:p4,即,Matlab 实现: function tout,yout=rk_4(odefile,tspan,y0) y0初值列向量t0=tspan(1); th=tspan(2);if length(tspan)=3, h=tspan(3); % tspan=t0,th,helse, h=tspan(2)-tspan(1); th=tspan(end); end 等间距数组 tout=t0:h:t
13、h; 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;yout=yout; y0;end 由效果看,该算法不是一个较好的方法。,5.2.3 一阶微分方程组的数值解 5.2.3.1 四阶五级Runge-Kutta-Felhberg算法,通
14、过误差向量调节步长,此为自动变步长方法。 四阶五级RKF算法有参量系数表。,5.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求解区间,x0初值问题的初始状态变量。,描述需要求解的微分方程组(Fun):方式1 不需附加变量function xd=funname(t,x)方式2 可以使用附加变量function xd=fun
15、name(t,x,flag,p1,p2,) t是时间变量或自变量(必须给);x为状态向量;xd为返回状态向量的导数;flag用来控制求解过程,指定初值,即使初值不用指定,也必须有该变量占位。,修改变量(options):方式1: odeset( ) 函数设置,例如:options=odeset(RelTol,1e-7);方式2:直接修改options成员变量,例如:options= odeset; options. RelTol= 1e-7;,例:,编写函数,function xdot = lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3);. -10*x(2)+10*
16、x(3) ;. -x(1)*x(2)+28*x(2)-x(3);,t_final=100; x0=0;0;1e-10; %输入计算终止时间和初值 t,x=ode45(lorenzeq,0,t_final,x0); plot(t,x),% 调用Matlab函数 figure; % 打开新图形窗口 plot3(x(:,1),x(:,2),x(:,3); axis(10 42 -20 20 -20 25); % 根据实际数值手动设置坐标系,可采用comet3( )函数绘制动画式的轨迹。 comet3(x(:,1),x(:,2),x(:,3),描述微分方程是常微分方程初值问题数值求解的关键。 f1=i
17、nline(-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);得出完全一致的结果。,5.2.3.3 MATLAB 下带有附加参数的微分方程求解,例:,编写函数 function xdot=lorenz1(t,x,flag,beta,rho,sigma)%
18、 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=3/8; r2=10; s2=28;%变量名不必一定为 等 t2,x2=ode45(lorenz1,0,t_final,x0,b2,r2,s2); plot(t2,x2), %options位置为,表示不需修改控制选项 figure; plot3(x2(:,1),x2(:,2),x2(:,3); axis(10 42 -20 20 -20 25);,例:求解,时的
19、Lorenz微分方程, 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); 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
20、,sigma); flag变量是不能省略,t_final=100; x0=0;0;1e-10; b2=2; r2=5; s2=20; t2,x2=ode45(f2,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); %得出的结果与前面一致,5.2.4 微分方程转换 5.2.4.1 单个高阶常微分方程处理方法,例:,函数描述为:function y=vdp_eq(t,x,flag,mu)y=x(2);
21、 -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);由于变步长所采用的步长过小,所需时间较长,导致输出的
22、y矩阵过大,超出计算机存储空间容量。所以不适合采用ode45()来求解,可用刚性方程求解算法ode15s( )。,5.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;
23、0; -1.04935751;% 输入初值 tic, t,y=ode45(apolloeq,0,20,x0); toc elapsed_time =0.5000 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.2970 length(t1), plot(y1(:,
24、1),y1(:,3), ans =1873, min(diff(t1) ans =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); toc elapsed_time =0.8590 plot(y1(:,1),y1(:,3) % 绘制出轨迹曲线显而易见,这样求解 是错误的,应该采用 更小的步长。, tic, t2,y2=rk_4(apolloeq,0,20,0.001,x0); toc elapsed_time =12.4380
25、 计算时间过长 plot(y2(:,1),y2(:,3) % 绘制出轨迹曲线严格说来某些点仍不 满足106的误差限, 所以求解常微分方程 组时建议采用变步长 算法,而不是定步长 算法。,例:试将二元方程组:,若两个高阶微分方程同时含有最高价导数项,需要先进行相应处理,再应用上述变换方法,先消去其中一个高阶导数,求解:,用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-5+x4*x1-x3)
26、/(3*x2+2*x4)dy =(2*x42*x1+5-x4*x1+x3)/(3*x2+2*x4)对于更复杂的问题来说,手工变换的难度将很大,所以如有可能,可采用计算机去求解有关方程,获得解析解。如不能得到解析解,也需要在描写一阶常微分方程组时列写出式子,得出问题的数值解。,5.3特殊微分方程的数值解 5.3.1 刚性微分方程的求解,刚性微分方程刚性过程:微分方程描述的变化过程中,若包含着多 个相互作用但变化速度相差十分悬殊的子过程,这样一类过程就认为具有 “刚性”。刚性微分方程:描述“刚性”过程的微分方程称为刚性微分方程,相应的初值问题称为“刚性问题” 。MATLAB采用求解函数ode15s
27、(),该函数的调用格式和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 =1.2970,作图 plot(t,y(:,1); figure; plot(t,y(:,2) Vandpol方程输出结果y=(y,y)y(:,1): 解曲线;部分曲线变化较平滑,某些点上变化在较快, y(
28、:,2) :解的变化率,例:定义函数 function dy=c5exstf2(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(c5exstf2,0,100,0;1); toc elapsed_time =33.0630 length(t2), plot(t2,y2) ans =356941,步长分析: format long, min(diff(t2), max(diff(t2) ans =0.00022220693884 0.0021497178718
29、4 plot(t2(1:end-1),diff(t2),方法二,用ode15s()代替ode45() opt=odeset; opt.RelTol=1e-6; tic,t1,y1=ode15s(c5exstf2,0,100,0;1,opt); toc elapsed_time =0.2340 length(t1), plot(t1,y1),legend(t1,y1) ans =169,5.3.2 隐式微分方程求解,隐式微分方程为不能转化为显式常微分方程组的方程 例:,编写函数: function dx=c5ximp(t,x)A=sin(x(1) cos(x(2); -cos(x(2) sin(
30、x(1);B=1-x(1); -x(2); dx=inv(A)*B; 求解: opt=odeset; opt.RelTol=1e-6; t,x=ode45(c5ximp,0,10,0; 0,opt); plot(t,x),应用matlab函数:ode15i,格式:,%如果不能直接确定初值,可用函数decic()得出相应初值,例:,5.3.3 微分代数方程求解,微分代数方程(differential algebra equation):指微分方程中某些变量间满足相应代数方程的约束,设,例:,编写函数function dx=c5eqdae(t,x)dx=-0.2*x(1)+x(2)*x(3)+0.
31、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(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(c5eqdae,0,20,x0,options); plot(t,x),编写函数: function dx=c5eqdae1(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*
32、x(2)*(1-x(1)-x(2)-2*x(2)*x(2);, x0=0.8; 0.1; fDae=inline(-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),t,x); t1,x1=ode45(fDae,0,20,x0); plot(t1,x1,t1,1-sum(x1),5.3.4延迟微分方程求解,sol:结构体数据,sol.x:时间向量t, sol.y:状态向量。,例:,编写函数: function dx=c5exdde(t,x,z)xlag1=z(:,1);
33、 %第一列表示提取xlag2=z(:,2);dx=1-3*x(1)-xlag1(2)-0.2*xlag2(1)3-xlag2(1);x(3); 4*x(1)-2*x(2)-3*x(3);,历史数据函数: function S=c5exhist(t)S=zeros(3,1);,求解: lags=1 0.5; tx=dde23(c5exdde,lags,zeros(3,1),0,10); plot(tx.x,tx.y(2,:) 与ode45()等返回的x矩阵不一样,它是按行排列的。,5.4边值问题的计算机求解,5.4.1 边值问题的打靶算法,将边值问题转化为初值问题考虑。或者说适当选择初始值使初值
34、问题的解满足边值条件。然后用求解初值问题的任一种有效的数值方法求解。,基本思路:,数学方法描述(以二阶方程为例),其相应边值条件,线性方程边值问题的打靶算法:,解法步骤,编写函数: 线性的 function t,y=shooting(f1,f2,tspan,x0f,varargin)t0=tspan(1); tfinal=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=(
35、gb-ga*y1(end,1)-yp(end,1)/y2(end,1);t,y=ode45(f2,tspan,ga;m,varargin);,例:编写函数: function xdot=c5fun1(t,x)xdot=x(2); -2*x(1)+3*x(2);function xdot=c5fun2(t,x)xdot=x(2); t-2*x(1)+3*x(2); t,y=shooting(c5fun1, c5fun2,0,1,1;2); plot(t,y),原方程的解析解为解的检验 y0=(exp(2)-3)*exp(t)+(3-exp(1)*exp(2*t)/(4*exp(1)*(exp(1
36、)-1)+3/4+t/2; norm(y(:,1)-y0) % 整个解函数检验 ans =4.4790e-008 norm(y(end,1)-2) % 终点条件检验 ans =2.2620e-008,非线性方程边值问题的打靶算法:,m:用Newton迭代法处理,编写函数: function t,y=nlbound(funcn,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,ts
37、pan,ga;m;0;1,varargin);m=m0-(v(end,1)-gb)/(v(end,3);endt,y=ode45(funcn,tspan,ga;m,varargin);,例:编写两个函数: function xdot=c5fun3(t,x)xdot=x(2); 2*x(1)*x(2); x(4); 2*x(2)*x(3)+2*x(1)*x(4);function xdot=c5fun4(t,x)xdot=x(2); 2*x(1)*x(2);, t,y=nlbound(c5fun4,c5fun3,0,pi/2,-1,1,1e-8); plot(t,y); set(gca,xlim
38、,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,5.4.2 线性微分方程的有限差分算法,把等式左边用差商表示。,编写函数: function x,y=fdiff(funcs,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
39、(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,i+1)=v(i); A(i+1,i)=w(i+1); endy=inv(A)*b; x=x tfinal; y=ga; y; gb;,例:编写函数: function y=c5fun5(x,key)switch keycase 1, y=1+x;case 2, y=1-x;otherwise, y=1+x.2;end t,y=fdiff(c5fun5,0,1,1,4,50); plot(t,y),