科學(xué)工程計算與matlab編程5

1、第五章 微分方程問題的解法 常系數(shù)線性微分方程的解析解方法 常微分方程問題的數(shù)值解法 微分方程問題算法概述 四階定步長 Runge-Kutta算法及 MATLAB 實現(xiàn) 一階微分方程組的數(shù)值解 微分方程轉(zhuǎn)換 特殊微分方程的數(shù)值解 邊值問題的計算機求解 偏微分方程的解5.1 常系數(shù)線性微分方程的解析解方法1212112( )( )( )( )( )( )nnnnnnnnd y tdy tdy tdy taaaa y tf tdtdtdtdt數(shù)學(xué)描述：特征方程（多項式代數(shù)方程）：121210nnnnnsa sa sasa該代數(shù)方程的根稱為原常系數(shù)方程的特征根根據(jù)特征根的情況可求得原方程的解析解 格

2、式： 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)*cos(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+

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

4、6*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) 如果用推導(dǎo)的方法求Ci的值，每個系數(shù)的解析解至少要寫出10數(shù)行,故可采用有理式近似 的方式表示. vpa(y,10) %有理近似值ans =1.196361839*e

5、xp(-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.9059*exp(-4.*t) 例：求解 x,y=dsolve(D2x+2*Dx=x+2*y-exp(-t), Dy=4*x+3*

6、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 line 292x =t-Int(1/(a-a3+1),a=.x)+C1=0故只有部分非線性微分方程可直接求

7、出其解析解。5.2 微分方程問題的數(shù)值解法算法概述 微分方程求解的誤差與步長問題：10000( ,)xxhf txR0( )xtfunction outx,outy=MyEuler(fun,x0,xt,y0,PointNum) % fun 表示f(x,y); x0,xt:自變量的初值和終值; y0:函數(shù)在x0處的值,其可以為向量形式; PointNum表示自變量在x0,xt上取的點數(shù)if nargin5 | PointNum=0 %PointNum 默認值為100 PointNum=100;endif narginf1=(x,y)sin(x)+y; 2*pi/0.4 % 計算 xt=2pi 時

8、應(yīng)取得點數(shù)ans = 15.7080 x1,y1=myeuler_1(f1,0,2*pi,1,16); %h=0.4 歐拉法所的的解x11,y11=MyEuler(myfun01,0,2*pi,1,32); %歐拉法所的的解h=0.2和h=0.4的數(shù)值解。y=dsolve(Dy=y+sin(t),y(0)=1); %該常微分方程的解析解for k=1:33t(k)=x11(k);y2(k)=subs(y,t(k); %求其對應(yīng)點的離散解endplot(x1,y1,+b,x11,y11,og,x11,y2,*r)legend(h=0.4的歐拉法解,h=0.2的歐拉解,符號解)誤差判斷：5.2 微

9、分方程問題的數(shù)值解法5.2.1 算法概述0.( , ), ( )dxf t xatbdtx ax 微分方程求解的誤差與步長問題：000121( )1,2,.,:nkkkkatx tattttbyknhtth時刻系統(tǒng)狀態(tài)向量的值為數(shù)值解法：即求微分方程初值問題的解在若干點 的近似值 （）的方法。步長（一般取等步長 ）應(yīng)用差商近似導(dǎo)數(shù)舍入誤差10000( ,)xxhf txR0( )x t提高數(shù)值解精度的方法之一： 減小步長h的值； function outx,outy=MyEuler(fun,x0,xt,y0,PointNum) % fun 表示f(x,y); x0,xt:自變量的初值和終值;

10、%y0:函數(shù)在x0處的值,其可以為向量形式; %PointNum表示自變量在x0,xt上取的點數(shù)if nargin=4 | PointNumf1=(x,y)sin(x)+y; 2*pi/0.4 % 計算 xt=2pi 時應(yīng)取得點數(shù)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的數(shù)值解。y=dsolve(Dy=y+sin(t),y(0)=1); %該常微分方程的解析解for k=1:33t(k)=x11(k)

11、;y2(k)=subs(y,t(k); %求其對應(yīng)點的離散解endplot(x1,y1,+b,x11,y11,og,x11,y2,*r)legend(h=0.4的歐拉法解,h=0.2的歐拉解,符號解)誤差判斷：改進的Euler算法 ( , )dxf t xdt1 kkttdt11()( )kkkkx tx txx1 kkttdt11( ,)(,)2kkkkhf txf tx1111( ,)( ,)(,)2kkkkkkkkkkxxhf txhxxf txf tx改進Euler迭代公式Euler公式梯形公式function Xout,Yout=MyEulerPro(fun,x0,xt,y0,Poi

12、ntNumber) %MyEulerPro 用改進的歐拉法解微分方程if nargin=4 | PointNumber5.2.2 四階定步長Runge-Kutta算法及Matlab實現(xiàn)Euler算法：一階Taylor多項式近似函數(shù)基本思想：2()O h誤差：RK一般算法：推廣：高級Taylor多項式近似函數(shù)提高數(shù)值解精度11111( ,)(,) (2,3,4,.)nniininijjpjnnjjkf txkf ta h xb kixxhkp為展開的項數(shù)例如：p4，即Matlab 實現(xiàn)：實現(xiàn)： function tout,yout=rk_4(odefile,tspan,y0) y0初值列向量 t

13、0=tspan(1); th=tspan(2); if length(tspan) comet3(x(:,1),x(:,2),x(:,3) 描述微分方程是常微分方程初值問題數(shù)值求解的關(guān)鍵。 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

14、25);得出完全一致的結(jié)果。5.2.3.3 MATLAB 下帶有附加參數(shù)的微分方程求解 例： 編寫函數(shù)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);3/8,10,28 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);

15、 plot(t2,x2), %options位置為，表示不需修改控制選項 figure; plot3(x2(:,1),x2(:,2),x2(:,3); axis(10 42 -20 20 -20 25);例：求解時的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(

16、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變量是不能省略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)

17、,x2(:,3); axis(0 72 -20 22 -35 40);%得出的結(jié)果與前面一致475.2.4 微分方程轉(zhuǎn)換5.2.4.1 單個高階常微分方程處理方法4849 例：1222121(1)xxxxxx 12;xyxy 函數(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,t_final,x0,mu); plo

18、t(t1,y1,t2,y2,:) figure; plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),:) 50 x0=2;0; t_final=3000; mu=1000; t,y=ode45(vdp_eq,0,t_final,x0,mu); 由于變步長所采用的步長過小，所需時間較長，導(dǎo)致輸出的y矩陣過大，超出計算機存儲空間容量。所以不適合采用ode45()來求解，可用剛性方程求解算法ode15s( )。 516.2.4.2 高階常微分方程組的變換方法52 例：53 54 描述函數(shù)： function dx=apolloeq(t,x) mu=1/82.45; mu1=1

19、-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;55 求解： x0=1.2; 0; 0; -1.04935751;% 輸入初值 tic, t,y=ode45(apolloeq,0,20,x0); tocelapsed_time = 0.5000 length(t), plot(y(:,1),y(:,3)ans = 689得出的

20、軌道不正確，默認精度RelTol設(shè)置得太大，從而導(dǎo)致的誤差傳遞，可減小該值。56 改變精度： options=odeset; options.RelTol=1e-6; tic, t1,y1=ode45(apolloeq,0,20,x0,options); tocelapsed_time = 0.2970 length(t1), plot(y1(:,1),y1(:,3),ans = 187357 min(diff(t1)ans = 1.8927e-004 plot(t1(1:end-1), diff(t1)58 例：59 x0=1.2; 0; 0; -1.04935751; tic, t1,y1

21、=rk_4(apolloeq,0,20,0.01,x0); tocelapsed_time = 0.8590 plot(y1(:,1),y1(:,3) % 繪制出軌跡曲線顯而易見，這樣求解是錯誤的，應(yīng)該采用更小的步長。 60 tic, t2,y2=rk_4(apolloeq,0,20,0.001,x0); tocelapsed_time = 12.4380 計算時間過長 plot(y2(:,1),y2(:,3) % 繪制出軌跡曲線嚴(yán)格說來某些點仍不滿足106的誤差限，所以求解常微分方程組時建議采用變步長算法，而不是定步長算法。 61 例：試將二元方程組：/2yyxx由第一個式子可得：若兩個高階

22、微分方程同時含有最高價導(dǎo)數(shù)項，需要先進行相應(yīng)處理，再應(yīng)用上述變換方法先消去其中一個高階導(dǎo)數(shù)求解：2102623xyxyxxyyx帶入第二個式子得：622314124422102623xxx xx x xxx所以： 1234,xx xxyxyx2431414425223xxx xx xxx 從而得： 63 用MATLAB符號工具箱求解，令1234,xx xx xy xy dxx dyy 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*

23、x1*x2-5+x4*x1-x3)/(3*x2+2*x4) dy = (2*x42*x1+5-x4*x1+x3)/(3*x2+2*x4) 對于更復(fù)雜的問題來說，手工變換的難度將很大，所以如有可能，可采用計算機去求解有關(guān)方程，獲得解析解。如不能得到解析解，也需要在描寫一階常微分方程組時列寫出式子，得出問題的數(shù)值解。645.3特殊微分方程的數(shù)值解 5.3.1 剛性微分方程的求解 剛性微分方程 剛性過程：微分方程描述的變化過程中，若包含著多 個相互作用但變化速度相差十分懸殊的子過程，這樣一類過程就認為具有 “剛性”。 剛性微分方程：描述“剛性”過程的微分方程稱為剛性微分方程，相應(yīng)的初值問題稱為“剛性

24、問題” 。 MATLAB采用求解函數(shù)ode15s()，該函數(shù)的調(diào)用格式和ode45()完全一致。t,x=ode15s(Fun,t0,tf,x0,options,p1,p2,) 65 例：計算 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 = 1.297066作圖 plot(t,y(:,1); figure; plot(t,y(:,2)Vandpol方程輸出結(jié)果y=(y,y) y(:,1):

25、 解曲線；部分曲線變化較平滑，某些點上變化在較快, y(:,2) ：解的變化率67 例：定義函數(shù)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;68 方法一 tic,t2,y2=ode45(c5exstf2,0,100,0;1); tocelapsed_time = 33.0630 length(t2), plot(t2,y2)ans = 35694169 步長分析： format long, min(diff(t2), max(diff(t2)ans

26、 = 0.00022220693884 0.00214971787184 plot(t2(1:end-1),diff(t2)70 方法二，用ode15s()代替ode45() opt=odeset; opt.RelTol=1e-6; tic,t1,y1=ode15s(c5exstf2,0,100,0;1,opt); tocelapsed_time = 0.2340 length(t1), plot(t1,y1), legend(t1,y1)ans = 169716.3.2 隱式微分方程求解 隱式微分方程為不能轉(zhuǎn)化為顯式常微分方程組的方程例：72 編寫函數(shù)：function dx=c5ximp(

27、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(c5ximp,0,10,0; 0,opt); plot(t,x),求解過程中Matlab若提示矩陣奇異的信息，則得出的數(shù)值解無意義.73應(yīng)用matlab函數(shù)：ode15i格式:%如果不能直接確定初值，可用函數(shù)decic()得出相應(yīng)初值74例：75765.3.3 微分代數(shù)方程求解微分代數(shù)方程(differential algebra equation)：指微分方程

28、中某些變量間滿足相應(yīng)代數(shù)方程的約束設(shè) ( , ):M( , ):f ttxx通常微分方程的非齊次項；奇異矩陣 為微分代數(shù)方程77例：78 編寫函數(shù) function dx=c5eqdae(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(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(c

29、5eqdae,0,20,x0,options); plot(t,x)注意這里使用ode45()函數(shù)求解會出錯79編寫函數(shù)：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*x(2)*(1-x(1)-x(2)-2*x(2)*x(2); 80 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);

30、t1,x1=ode45(fDae,0,20,x0); plot(t1,x1,t1,1-sum(x1)注意這里使用ode45()函數(shù)求解不會出錯也可應(yīng)用inline()函數(shù)定義815.3.3 微分代數(shù)方程求解 微分代數(shù)方程(differential algebra equation)：指微分方程中某些變量間滿足相應(yīng)代數(shù)方程的約束 設(shè) ( , ):M( , ):massf ttxx通常微分方程的非齊次項；奇異矩陣,即為質(zhì)量矩陣（） 為微分代數(shù)方程82例：83 編寫函數(shù) function dx=c5eqdae(t,x) dx=-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2); 2*

31、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微分代數(shù)方程中的質(zhì)量矩陣（控制參數(shù)） x0=0.8; 0.1; 0.1; t,x=ode15s(c5eqdae,0,20,x0,options); plot(t,x)注意這里使用ode45()函數(shù)求解會會出錯84編寫函數(shù)：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

32、)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2); 85 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)注意這里使用ode45()函數(shù)求解不會不會出錯也可應(yīng)用inline()函數(shù)定義865.3.4延遲微分方程求解sol:結(jié)構(gòu)體數(shù)據(jù)，sol.x:時間向量t, sol.y:狀態(tài)向量。

33、 1120:,:nfftt延 遲 微 分 方 程 ，時 的 狀 態(tài) 變 量 值 函 數(shù) 。87 例： 88編寫函數(shù)：function dx=c5exdde(t,x,z) xlag1=z(:,1); %第一列表示提取 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);1( )x歷史數(shù)據(jù)函數(shù)：function S=c5exhist(t) S=zeros(3,1);89 求解： lags=1 0.5; tx=dde23(c5exdde,lags,zeros(3,1),0,10);

34、 plot(tx.x,tx.y(2,:)與ode45()等返回的x矩陣不一樣，它是按行排列的。905.4邊值問題的計算機求解915.4.1 邊值問題的打靶算法將邊值問題轉(zhuǎn)化為初值問題考慮?；蛘哒f適定選將邊值問題轉(zhuǎn)化為初值問題考慮?；蛘哒f適定選擇初始值使初擇初始值使初值問題的解滿足邊值條件。然后用值問題的解滿足邊值條件。然后用求解初值問題的任一種有效的數(shù)值方法求解。求解初值問題的任一種有效的數(shù)值方法求解。 基本思路：基本思路：( , ,)( , ,)( ), ( )( ),( )yF x y yyF x y yy ay by ay am1m2m123m92數(shù)學(xué)方法描述（以二階方程為例）( ),

35、( )aby ay b其相應(yīng)邊值條件線性方程邊值問題的打靶算法：93( )( ) ( )( ) ( )( )( ),( )y xp x y xq x y xf xy ay am初值問題解法步驟11).( )( ) ( )( ) ( )0, ( )1, ( )0 ( )y xp x y xq x y xy ay ay b求 的數(shù)值解22).( )( ) ( )( ) ( )0, ( )0, ( )1 ( )y xp x y xq x y xy ay ay b求 的數(shù)值解33).( )( ) ( )( ) ( )( ), ( )0, ( )0 ( )y xp x y xq x y xf xy ay

36、 ay b求 的數(shù)值解13224).bam 若0，求 ( ), ( )aby ay b5).邊值問題94編寫函數(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(en

37、d,1); t,y=ode45(f2,tspan,ga;m,varargin);( , ,)F x y y12;xyxy95例：編寫函數(shù)：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)1212221;32xy xyxxxxx96原方程的解析解為解的檢驗 y0=(exp(2)-3)*exp(t)+(3-exp(1)*exp(2*t)/(4*

38、exp(1)*(exp(1)-1)+3/4+t/2; norm(y(:,1)-y0) % 整個解函數(shù)檢驗ans = 4.4790e-008 norm(y(end,1)-2) % 終點條件檢驗ans = 2.2620e-00897 非線性方程邊值問題的打靶算法：113(), ()biivbmmvb( , ,)( , ,)( ), ( )( ),( )abayF x y yyF x y yy ay by ay amm:用Newton迭代法處理1324(; ),(; ),(; ),(; )iiiiyvy m xvm xmyvy m xvm xm98編寫函數(shù)：function t,y=nlbound(

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

40、 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);1221234121243412( ,)( ,)( ,)vvvF x v vvvF x v vF x v vvvvvv100 t,y=nlbound(c5fun4,c5fun3,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

41、.2815e-0061016.4.2 線性微分方程的有限差分算法把等式左邊用差商表示。( ), ( )aby ay b( , ,)yF x y y112iiiyyyh102103編寫函數(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 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+

42、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;104例：編寫函數(shù)：function y=c5fun5(x,key) switch key case 1, y=1+x; case 2, y=1-x; otherwise, y=1+x.2; end t,y=fdiff

43、(c5fun5,0,1,1,4,50); plot(t,y)1056.5 偏微分方程求解入門 6.5.1 偏微分方程組求解函數(shù)描述：106 邊界條件的函數(shù)描述： 初值條件的函數(shù)描述： u0=pdeic(x)107 例：108 函數(shù)描述： 109function 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*du(1); 0.17*du(2); 110描述邊界條件的函數(shù)function pa,qa,pb,qb=c7mpbc(xa,ua,xb,ub,t) pa=0

44、; ua(2); qa=1;0; pb=ub(1)-1; 0; qb=0;1;120.024/0.17/uxfux111 描述初值：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)1126.5.2 二階偏微分方程的數(shù)學(xué)描述 橢圓型偏微分方程：113 拋物線型偏微分方程： 雙曲型偏微分方程： 特征值型偏微分方程：1146.5.3 偏微分方程的求解界面

45、應(yīng)用簡介 6.5.3.1 偏微分方程求解程序概述 啟動偏微分方程求解界面 在 MATLAB 下鍵入 pdetool 該界面分為四個部分 菜單系統(tǒng) 工具欄 集合編輯 求解區(qū)域1156.5.3.2 偏微分方程求解區(qū)域繪制1）用工具欄中的橢圓、矩形等繪制一些區(qū)域。2）在集合編輯欄中修改其內(nèi)容。 如（R1E1E2）E33）單擊工具欄中 按紐可得求解邊界。4）選擇Boundary-Remove All Subdomain Borders菜單項，消除相鄰區(qū)域中間的分隔線。5）單擊 按紐可將求解區(qū)域用三角形劃分成網(wǎng)格?？捎?按紐加密。1166.5.3.3 偏微分方程邊界條件描述選擇Boundary-Spec

46、ify Boundary Conditions菜單狄利克雷條件，諾伊曼條件。1176.5.3.4 偏微分方程求解舉例 例：求解：1）繪制求解區(qū)域。2）描述邊界條件（Boundary-Specify Boundary Conditions）。3）選擇偏微分方程的類型。單擊工具欄中的PDE圖標(biāo)，在打開的新窗口選擇Hyperbolic選項，輸入?yún)?shù)c,a,f,d.4）求解。單擊工具欄中的等號按鈕。 118顯示：1）圖形顏色表示t=0時u(x,y)的函數(shù)值。2）單擊工具欄中的三維圖標(biāo)將打開一新的對話框，若再選擇Contour可繪制等值線，若選擇Arrows選項可繪制引力線。若單獨選擇Height（3d

47、-plot)，則在另一窗口繪制出三維圖形。3）可在單擊三維圖標(biāo)打開的新對話框中，對Property欄目的各個項目重新選擇。4）可修改微分方程的邊界條件，重新求解。動畫：1）Solve-Parameters對話框時間向量改為0:0.1:2。2）三維圖標(biāo)打開的對話框中選擇Animation選項，單擊Options按紐設(shè)置播放速度。Plot-Export Movie 菜單。1196.5.3.5 函數(shù)參數(shù)的偏微分方程求解 例：(橢圓型）120 求解：1）求解區(qū)域不變。2）描述邊界條件，u=0。3）選擇偏微分方程的類型。單擊工具欄中的PDE圖標(biāo)，在打開的新窗口選擇Elliptic選項，輸入?yún)?shù)c=1./

48、sqrt(1+ux.2+uy.2), a=x.2+y.2 , f=exp(-x.2-y.2).4)再打開Solve-Parameters對話框，選定Use nonlinear solve屬性（該屬性只適于橢圓性偏微分方程）5）求解。單擊工具欄中的等號按鈕。121例：橢圓微分方程的拉普拉斯形式：其自變量取值范圍D為： 邊界條件為：22222( , )( , )( , )0u x yu x yu x yxy04,04xy22( ,0),( ,4)16cos( ),(0, ),(4, )16cos( ).u xxu xxuyyuyy1221、調(diào)整坐標(biāo)范圍:options-Axes Limits2、設(shè)

49、定矩形區(qū)域，雙擊確定調(diào)整。3、設(shè)置邊界條件：單擊邊界按鈕 ，并雙擊相應(yīng)的邊彈出對話框，設(shè)定邊界條件。4、設(shè)定方程。按“PDE”按鈕選定Elliptic.5、單擊 按紐可將求解區(qū)域用三角形劃分成網(wǎng)格。可用 按紐加密。6、單擊三維圖標(biāo)，設(shè)置所畫曲線特性。7、單擊“plot”或工具欄中的“=”按鈕。5.5 偏微分方程求解入門5.5.1 偏微分方程組求解1+1維初-邊值問題的數(shù)值解: pdepe( )pdepe( )可求解如下形式的方程0, , 0, 1, 2ftttaxb mor 初值條件（對任意x成立）：邊界條件： 其中：u0=funic(u0=funic(x) ),funbc(, )ababaa

50、bbpp q qx ux u t01, 0 xt例：0 函數(shù)描述function c,f,s=c5mpde(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*du(1); 0.17*du(2); 描述邊界條件的函數(shù)function pa,qa,pb,qb=c5mpbc(xa,ua,xb,ub,t) pa=0; ua(2); qa=1;0; pb=ub(1)-1; 0; qb=0;1;120.024/0.17/uxfux 描述初值：function u0=c5mpic(x) u0=1; 0;求解： x=0:.05:1; t=0:0.05:2; m=0; sol=pdepe(m,c5mpde,c5mpic,c5mpbc,x,t); surf(x,t,sol(:,:,1)， figure; surf(x,t,sol(:,:,2)05.5.2 二階偏微分方程的數(shù)學(xué)描述 橢圓型偏微分方程： 拋物線型偏微分方程：雙曲型偏微分方程：當(dāng)c為常數(shù)時：當(dāng)c為常數(shù)時： 特征值