MALAB微分方程問題的解法ppt課件

1、第六章 微分方程問題的解法 微分方程的解析解方法 常微分方程問題的數(shù)值解法 微分方程問題算法概述 四階定步長 Runge-Kutta算法及 MATLAB 實現(xiàn) 一階微分方程組的數(shù)值解 微分方程轉(zhuǎn)換 特殊微分方程的數(shù)值解 邊值問題的計算機求解 偏微分方程的解6.1 微分方程的解析解方法 格式： y=dsolve(f1, f2, , fm) 格式：指明自變量 y=dsolve(f1, f2, , fm ,x) fi即可以描述微分方程，又可描述初始條件或邊界條件。如： 描述微分方程時 描述條件時 (4)( )747ytDy(2)32 (2)3yD y例： syms t; u=exp(-5*t)*co

2、s(2*t+1)+5; uu=5*diff(u,t,2)+4*diff(u,t)+2*uuu =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(

3、0)=3, Dy(0)=2, D2y(0)=0, D3y(0)=0)分別處理系數(shù)，如： n,d=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)ans = -8704 185 % rat()最接近有理數(shù)的分數(shù)判斷誤差： 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)*s

4、in(2*t+1) + . 10,y(0)=1/2,Dy(pi)=1,D2y(2*pi)=0,Dy(2*pi)=1/5); 如果用推導(dǎo)的方法求Ci的值，每個系數(shù)的解析解至少要寫出10數(shù)行,故可采用有理式近似 的方式表示. 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.

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=dsolve(Dx=x*(1-x2)+1) Warning: Explicit solut

6、ion 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 故只有部分非線性微分方程有解析解。6.2 微分方程問題的數(shù)值解法6.2.1 微分方程問題算法概述微分方程求解的誤差與步長問題：6.2.2 四階定步長Runge-Kutta算法 及 MATLAB 實現(xiàn) function tout,yout=rk_4(odefile,tspan,y0) y0初值列向量 t0=tspan(1); th=

7、tspan(2); if length(tspan) t_final=100; x0=0;0;1e-10; % t_final為設(shè)定的仿真終止時間 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); % 根據(jù)實際數(shù)值手動設(shè)置坐標系 可采用comet3( )函數(shù)繪制動畫式的軌跡。 comet3(x(:,1),x(:,2),x(:,3) 描述微分方程是常微分方程初值問題數(shù)值求解的關(guān)鍵。 f1=inline(-8/3*

8、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); 得出完全一致的結(jié)果。6.2.3.3 MATLAB 下帶有附加參數(shù)的微分方程求解 例： 編寫函數(shù) function xdot=lorenz1(t,x,flag,beta,rho,sigma) flag變量是不能省

9、略的 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); axis(0 72 -20 22 -35 40); f2=inline(-beta*x(1)+x(2)

10、*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 微分方程轉(zhuǎn)換6.2.4.1 單個高階常微分方程處理方法 例： 函數(shù)描述為： 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(vdp_eq,0,t_final,x0,mu); mu=2; t2,y2=ode45(vdp_eq,0,

11、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); 由于變步長所采用的步長過小，所需時間較長，導(dǎo)致輸出的y矩陣過大，超出計算機存儲空間容量。所以不適合采用ode45()來求解，可用剛性方程求解算法ode15s( )。6.2.4.2 高階常微分方程組的變換方法 例： 描述函數(shù)： function dx=apolloeq(t,x) mu=1/82.45;

12、 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.8310 length(t), plot(y(:,1),y(:,3) ans = 689 得出

13、的軌道不正確， 默認精度RelTol設(shè)置 得太大，從而導(dǎo)致的 誤差傳遞，可減小該 值。 改變精度： 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) 例： x0=1.2; 0; 0; -1.04935751; tic, t1,y1

14、=rk_4(apolloeq,0,20,0.01,x0); toc elapsed_time = 4.2570 plot(y1(:,1),y1(:,3) % 繪制出軌跡曲線 顯而易見，這樣求解 是錯誤的，應(yīng)該采用 更小的步長。 tic, t2,y2=rk_4(apolloeq,0,20,0.001,x0); toc elapsed_time = 124.4990 計算時間過長 plot(y2(:,1),y2(:,3) % 繪制出軌跡曲線 嚴格說來某些點仍不 滿足106的誤差限， 所以求解常微分方程 組時建議采用變步長 算法，而不是定步長 算法。 例：用MATLAB符號工具箱求解，令 syms

15、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) 對于更復(fù)雜的問題來說，手工變換的難度將很大，所以如有可能，可采用計算機去求解有關(guān)方程，獲得解析解。如不能得到解析解，也需要在描寫一階常微分方程組時列寫出式子，得出問題的數(shù)值解。1234,xx xx xy xy dxx dyy6.3特殊微分方程的數(shù)值解 6.3.1

16、剛性微分方程的求解 剛性微分方程 一類特殊的常微分方程，其中一些解變化緩慢，另一些變化快，且相差懸殊，這類方程常常稱為剛性方程。 MATLAB采用求解函數(shù)ode15s()，該函數(shù)的調(diào)用格式和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.5240作圖 plot(t,

17、y(:,1); figure; plot(t,y(:,2) y(:,1)曲線變化較平滑, y(:,2)變化在某些點上較快。 例： 定義函數(shù) 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, min(diff

18、(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 = 1696.3.2 隱式微分方程求解 隱式微分方程為不能轉(zhuǎn)化為顯式常微分方程組的方程 例： 編寫函數(shù)： func

19、tion 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 微分代數(shù)方程求解例： 編寫函數(shù) 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(

20、3)-1; M=1,0,0; 0,1,0; 0,0,0; options=odeset; options.Mass=M; Mass微分代數(shù)方程中 的質(zhì)量矩陣控制參數(shù)） x0=0.8; 0.1; 0.1; t,x=ode15s(c7eqdae,0,20,x0,options); plot(t,x)編寫函數(shù)：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.2*

21、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:結(jié)構(gòu)體數(shù)據(jù)，sol.x:時間向量t, sol.y:狀態(tài)向量。 1120:,:nfftt延 遲 微 分 方 程 ，時 的 狀 態(tài) 變 量 值 函 數(shù) 。 例：編寫函數(shù)：function dx=c7exdde(t,x,z) xlag1=z(:,1); %第一列表示提取 xlag2=z(:,2

22、); 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);歷史數(shù)據(jù)函數(shù)：function S=c7exhist(t) S=zeros(3,1);1( )x 求解： 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 邊值問題的打靶算法數(shù)學方法描述： 以二階方程為例 ( , ,)( , ,)( ), ( )( ),( )yF

23、 x y yyF x y yy ay by ay am2131121()mmmm1m2m12m編寫函數(shù)： 線性的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=(gb-ga*y1(end,1)-yp(end,1)/y2(end,1);

24、 t,y=ode45(f2,tspan,ga;m,varargin);( , ,)F x y y例：編寫函數(shù)：function xdot=c7fun1(t,x) xdot=x(2); -2*x(1)+3*x(2);function xdot=c7fun2(t,x) xdot=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(:,

25、1)-y0) % 整個解函數(shù)檢驗ans = 4.4790e-008 norm(y(end,1)-2) % 終點條件檢驗ans = 2.2620e-008非線性方程邊值問題的打靶算法： 用Newton迭代法處理(,)yFxyy( , ,)( , ,)( ), ( )( ),( )yF x y yyF x y yy ay by ay am編寫函數(shù)：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 (n

26、orm(m-m0)tol), m0=m; t,v=ode45(funcv,tspan,ga;m;0;1,varargin); m=m0-(v(end,1)-gb)/(v(end,3); end t,y=ode45(funcn,tspan,ga;m,varargin);例：編寫兩個函數(shù)：function xdot=c7fun3(t,x) xdot=x(2); 2*x(1)*x(2); x(4); 2*x(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

27、,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-0066.4.2 線性微分方程的有限差分算法把等式左邊用差商表示。( ), ( )aby ay b編寫函數(shù)：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

28、 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,i+1)=v(i); A(i+1,i)=w(i+1); end y=inv(A)*b; x=x tfinal; y=ga; y; gb;例：編寫函數(shù)：funct

29、ion y=c7fun5(x,key) switch key case 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 偏微分方程組求解函數(shù)描述： 邊界條件的函數(shù)描述： 初值條件的函數(shù)描述： u0=pdeic(x) 例： 函數(shù)描述： function c,f,s=c7mpde(x,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

30、*du(1); 0.17*du(2); 描述邊界條件的函數(shù) 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); surf(x,t,sol(:,:,1)， figure; surf(x,t,sol(:,:,2)6.5.2 二階偏微分方程的數(shù)學描述 橢圓型偏微分方程： 拋物線型偏微分方程： 雙曲型偏微分方程： 特征值型偏微分