微分方程問題的解法

1、第六章 微分方程問題的解法,微分方程的解析解方法常微分方程問題的數(shù)值解法微分方程問題算法概述四階定步長(zhǎng) Runge-Kutta算法及 MATLAB 實(shí)現(xiàn)一階微分方程組的數(shù)值解微分方程轉(zhuǎn)換特殊微分方程的數(shù)值解邊值問題的計(jì)算機(jī)求解偏微分方程的解,6.1 微分方程的解析解方法,格式： y=dsolve(f1, f2, …, fm)格式：指明自變量 y=dsolve(f

2、1, f2, …, fm ,’x’) fi即可以描述微分方程，又可描述初始條件或邊界條件。如： 描述微分方程時(shí) 描述條件時(shí),例：>> 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

3、)*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=&

4、#39;,... '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*

5、sin(1)-69/2)))]ans = -8704 185 % rat()最接近有理數(shù)的分?jǐn)?shù)判斷誤差：>> vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185)ans =.114731975864790922564144636e-4,>> y=dsolve(['D4y+10*D3y+35*D

6、2y+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的值，每個(gè)系數(shù)的解析解至少要寫

7、出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)

8、*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

9、x>> x=dsolve('Dx=x*(1-x^2)')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-x^2)+1')Warning: Explicit solution could not be found; implicit s

10、olution returned.> In D:\MATLAB6p5\toolbox\symbolic\dsolve.m at line 292x =t-Int(1/(a-a^3+1),a=``..x)+C1=0故只有部分非線性微分方程有解析解。,6.2 微分方程問題的數(shù)值解法6.2.1 微分方程問題算法概述,微分方程求解的誤差與步長(zhǎng)問題：,function [outx,outy]=MyEuler(fun,x0,xt,y

11、0,PointNum) % fun 表示f(x,y); x0,xt:自變量的初值和終值; y0:函數(shù)在x0處的值,其可以為向量形式; PointNum表示自變量在[x0,xt]上取的點(diǎn)數(shù)if nargin<5 | PointNum<=0 %PointNum 默認(rèn)值為100 PointNum=100;endif nargin<4

12、 % y0默認(rèn)值為0 y0=0;end,h=(xt-x0)/PointNum; %計(jì)算步長(zhǎng)hx=x0+[0:PointNum]'*h; %自變量數(shù)組y(1,:) = y0(:)'; %將輸入存為行向量，輸入為列向量形式for k = 1:PointNum f=feval(fun,x(k),y(k,:)); f=f(:)'; %計(jì)算f(x,y

13、)在每個(gè)迭代點(diǎn)的值 y(k + 1,:) =y(k,:) +h*f; %對(duì)于所取的點(diǎn)x迭代計(jì)算y值endouty=y; outx=x; plot(x,y) %畫出方程解的函數(shù)圖,例：求 h=0.2和h=0.4的數(shù)值解。function ex0801()[x1,y1]=MyEuler('myfun01',0

14、,2*pi,1,16); %歐拉法所的的解h1=2*pi/15 %計(jì)算取16點(diǎn)的步長(zhǎng)[x11,y11]=MyEuler('myfun01',0,2*pi,1,32); %歐拉法所的的解h2=pi/15 %計(jì)算取32點(diǎn)的步長(zhǎng),y=dsolve('Dy=y+sin(t)','y(0)=1');

15、 %該常微分方程的解析解for k=1:33t(k)=x11(k);y2(k)=subs(y,t(k)); %求其對(duì)應(yīng)點(diǎn)的離散解endplot(x1,y1,'+b',x11,y11,'og',x11,y2,'*r')legend('h=0.4的歐拉法解','h=0.2的歐拉解','符號(hào)解&#

16、39;),function f=myfun01(x,y) f=sin(x)+y;>> ex0801h1 = 0.4189h2 = 0.2094,function [Xout,Yout]=MyEulerPro(fun,x0,xt,y0,PointNumber) %MyEulerPro 用改進(jìn)的歐拉法解微分方程if nargin<5 | PointNu

17、mber<=0 %PointNumer默認(rèn)值為100 PointNumer=100;endif nargin<4 %y0默認(rèn)值為0 y0=0;end,h=(xt-x0)/PointNumber; %計(jì)算所取的兩離散點(diǎn)之間的距離x=x0+[0:PointNumber]'*h; %表示出離散的自變量xy(1,:)=y0

18、(:)';for i=1:PointNumber %迭代計(jì)算過程 f1=h*feval(fun,x(i),y(i,:)); f1=f1(:)'; f2=h*feval(fun,x(i+1),y(i,:)+f1); f2=f2(:)'; y(i+1,:)=y(i,:)+1/2*(f1+f2);endXout=x;Yout=y;,function ex0802()%比

19、較改進(jìn)歐拉法，簡(jiǎn)單歐拉方法以及微分方程符號(hào)解[x3,y3]=MyEulerPro('myfun01',0,2*pi,1,128);[x,y1]=MyEuler('myfun01',0,2*pi,1,128);%歐拉法所的的解y=dsolve('Dy=y+sin(t)','y(0)=1');%該常微分方程的符號(hào)解for k=1:129t(k)=x(k);y2(k)

20、=subs(y,t(k));%求其對(duì)應(yīng)點(diǎn)的離散解endplot(x,y1,'-b',x3,y3,'og',x,y2,'*r')legend('簡(jiǎn)單歐拉法解','改進(jìn)歐拉解','符號(hào)解'),>> ex0802,6.2.2 四階定步長(zhǎng)Runge-Kutta算法 及 MATLAB 實(shí)現(xiàn),function [tout,y

21、out]=rk_4(odefile,tspan,y0) ％y0初值列向量 t0=tspan(1); th=tspan(2); if length(tspan)<=3, h=tspan(3); % tspan=[t0,th,h] else, h=tspan(2)-tspan(1); th=tspan(end); end ％等間距數(shù)組 tout=[t0:h:th]'; yout=[]; f

22、or t=tout' k1=h*eval([odefile ‘(t,y0)’]); % odefile是一個(gè)字符串變量，為表示微分方程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([od

23、efile '(t+h,y0+k3)']); y0=y0+(k1+2*k2+2*k3+k4)/6; yout=[yout; y0']; end ％由效果看，該算法不是一個(gè)較好的方法。,6.2.3 一階微分方程組的數(shù)值解6.2.3.1 四階五級(jí)Runge-Kutta-Felhberg算法,通過誤差向量調(diào)節(jié)步長(zhǎng)，此為自動(dòng)變步長(zhǎng)方法。四階五級(jí)RKF算法有參量系數(shù)表。,

24、6.2.3.2 基于 MATLAB 的微分方程,求解函數(shù)格式1： 直接求解 [t,x]=ode45(Fun,[t0,tf],x0) 格式2： 帶有控制參數(shù) [t,x]=ode45(Fun,[t0,tf],x0，options)格式3： 帶有附加參數(shù)[t,x]=ode45(Fun,[t0,tf],x0,options,p1,p2,…)[t0,tf]求解區(qū)間，x0初值問題的初始狀態(tài)變量。,描述需要

25、求解的微分方程組： 不需附加變量的格式 function xd=funname(t,x) 可以使用附加變量 function xd=funname(t,x,flag,p1,p2,…) ％ t是時(shí)間變量或自變量（必須給），x為狀態(tài)向量， xd為返回狀態(tài)向量的導(dǎo)數(shù)。flag用來控制求解過程,指定初值，即使初值不用指定，也必須有該變量占位。修改變量：options唯一結(jié)構(gòu)體變量，用odes

26、et( )修改。 options=odeset(‘RelTol’,1e-7); options= odeset; options. RelTol= 1e-7;,例：自變函數(shù) function xdot = 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)];,&g

27、t;> t_final=100; x0=[0;0;1e-10]; % t_final為設(shè)定的仿真終止時(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

28、 -20 25]); % 根據(jù)實(shí)際數(shù)值手動(dòng)設(shè)置坐標(biāo)系,可采用comet3( )函數(shù)繪制動(dòng)畫式的軌跡。>> 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

29、)]'],'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é)果

30、。,6.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);

31、 -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位置為[]，

32、表示不需修改控制選項(xiàng)>> 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','

33、;x','flag','beta','rho','sigma');％ flag變量是不能省略的,6.2.4 微分方程轉(zhuǎn)換6.2.4.1 單個(gè)高階常微分方程處理方法,例：函數(shù)描述為： function y=vdp_eq(t,x,flag,mu) y=[x(2); -mu*(x(1)^2-1)*x(2)-x(1)];>> x0=[-0

34、.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; plo

35、t(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); 由于變步長(zhǎng)所采用的步長(zhǎng)過小，所需時(shí)間較長(zhǎng)，導(dǎo)致輸出的y矩陣過大，超出計(jì)算機(jī)存儲(chǔ)空間容量。所以不適合采用ode45()來求解，可用剛

36、性方程求解算法ode15s( )。,6.2.4.2 高階常微分方程組的變換方法,例：,描述函數(shù)： 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)-m

37、u1*(x(1)+mu)/r1^3-mu*(x(1)-mu1)/r2^3; x(4); -2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];,求解：>> x0=[1.2; 0; 0; -1.04935751];>> tic, [t,y]=ode45('apolloeq',[0,20],x0); tocelapsed_ti

38、me = 0.8310>> length(t), >> plot(y(:,1),y(:,3))ans = 689得出的軌道不正確，默認(rèn)精度RelTol設(shè)置得太大，從而導(dǎo)致的誤差傳遞，可減小該值。,改變精度：>> options=odeset; options.RelTol=1e-6;>> tic, [t1,y1]=ode45('apolloeq&#

39、39;,[0,20],x0,options); tocelapsed_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)),例：,&

40、gt;> 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)) % 繪制出軌跡曲線顯而易見，這樣求解是錯(cuò)誤的，應(yīng)該采用更小的步長(zhǎng)。,>> tic, [

41、t2,y2]=rk_4('apolloeq',[0,20,0.001],x0); tocelapsed_time = 124.4990 ％計(jì)算時(shí)間過長(zhǎng)>> plot(y2(:,1),y2(:,3)) % 繪制出軌跡曲線嚴(yán)格說來某些點(diǎn)仍不滿足10－6的誤差限，所以求解常微分方程組時(shí)建議采用變步長(zhǎng)算法，而不是定步長(zhǎng)算法。,例：,用MATLAB符號(hào)工具箱求解

42、，令 ％>> 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-

43、5)/(2*x4+3*x2)dy = (2*x4^2*x1-x4*x1+x3+5)/(2*x4+3*x2) 對(duì)于更復(fù)雜的問題來說，手工變換的難度將很大，所以如有可能，可采用計(jì)算機(jī)去求解有關(guān)方程，獲得解析解。如不能得到解析解，也需要在描寫一階常微分方程組時(shí)列寫出式子，得出問題的數(shù)值解。,6.3特殊微分方程的數(shù)值解 6.3.1 剛性微分方程的求解,剛性微分方程 一類特殊的常微分方程，其中一

44、些解變化緩慢，另一些變化快，且相差懸殊，這類方程常常稱為剛性方程。 MATLAB采用求解函數(shù)ode15s()，該函數(shù)的調(diào)用格式和ode45()完全一致。[t,x]=ode15s(Fun,[t0,tf],x0,options,p1,p2,…),例：％計(jì)算>> h_opt=odeset; h_opt.RelTol=1e-6; >> x0=[2;0]; t_final=3000;>>

45、; tic, mu=1000; [t,y]=ode15s('vdp_eq',[0,t_final],x0,h_opt,mu); tocelapsed_time = 2.5240,％作圖>> plot(t,y(:,1)); figure; plot(t,y(:,2)) y(:,1)曲線變化較平滑, y(:,2)變化在某些點(diǎn)上較快。,例：定義函數(shù)function dy=c7e

46、xstf2(t,y) dy=[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2; -10^4*y(1)+3000*(1-y(2))^2];,方法一>> tic,[t2,y2]=ode45('c7exstf2',[0,100],[0;1]); tocelapsed_time = 229.4700>> length(t2), pl

47、ot(t2,y2)ans = 356941,步長(zhǎng)分析：>> format long, [min(diff(t2)), max(diff(t2))]ans = 0.00022220693884 0.00214971787184>> plot(t2(1:end-1),diff(t2)),方法二，用ode15s()代替ode45()>> opt=odeset; opt.Rel

48、Tol=1e-6;>> tic,[t1,y1]=ode15s('c7exstf2',[0,100],[0;1],opt); tocelapsed_time = 0.49100000000000>> length(t1),>> plot(t1,y1)ans = 169,6.3.2 隱式微分方程求解,隱式微分方程為不能轉(zhuǎn)化為顯式常微分方程組的方程例：,編寫函數(shù)：f

49、unction 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.

50、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(3)-1];>> M=[1,0,0; 0,1,0; 0,0,0]; >> options=odeset;

51、 >> 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)

52、*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*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)]

53、9;],'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:時(shí)間向量t, sol.y:狀態(tài)向量。,例：,編寫函數(shù)???：function dx=c7exdde(t,x,z) xlag1=z(:,1); %第一列表示提取 xl

54、ag2=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)];歷史數(shù)據(jù)函數(shù)：function S=c7exhist(t) S=zeros(3,1);,求解：>> lags=[1 0.5]; tx=dde23('c7exdde',lags,zeros(3,1),[

55、0,10]);>> plot(tx.x,tx.y(2,:))％與ode45()等返回的x矩陣不一樣，它是按行排列的。,6.4邊值問題的計(jì)算機(jī)求解,6.4.1 邊值問題的打靶算法,數(shù)學(xué)方法描述： 以二階方程為例,,,,,,,,編寫函數(shù)： 線性的function [t,y]=shooting(f1,f2,tspan,x0f,varargin) t0=tspan(1); tfinal=

56、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,tsp

57、an,[ga;m],varargin);,例：編寫函數(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],[

58、1;2]); plot(t,y),原方程的解析解為解的檢驗(yàn)>> 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) % 整個(gè)解函數(shù)檢驗(yàn)ans = 4.4790e-008>> norm(y(end,1)-2) % 終點(diǎn)條件檢驗(yàn)ans = 2.26

59、20e-008,非線性方程邊值問題的打靶算法： 用Newton迭代法處理 m,編寫函數(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 (norm(m-m0)>tol), m0=m;

60、[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);,例：編寫兩個(gè)函數(shù)：function xdot=c7fun3(t,x) xdot=[x(2); 2*x(1)*x(2); x(4); 2*x(2)*x(3)

61、+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]);精確解：檢驗(yàn)：>&

62、gt; 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 線性微分方程的有限差分算法,把等式左邊用差商表示。,,編寫函數(shù)：function [x,y]=fdiff(funcs,tspan,x0f,n) t0=tspan(1);tfinal=tspan(2)

63、; 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+h^2*feval(funcs,x0,2); tmp=feval(funcs,x0,1); v=1+h*tmp/2; w=1-h*tmp/2; b=h^2*feval(funcs,x0,3); b(1)=b(1)-w(1)*g

64、a; 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ù)：function y=c7fun5(x,key) switch key case 1, y=1+x;

65、 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ù)描述：

66、 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*du(1); 0.17*du(2)];,描述邊界條件的函數(shù)function [pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t)

67、 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))， fi

68、gure; surf(x,t,sol(:,:,2)),6.5.2 二階偏微分方程的數(shù)學(xué)描述,橢圓型偏微分方程：,拋物線型偏微分方程：雙曲型偏微分方程：特征值型偏微分方程：,6.5.3 偏微分方程的求解界面應(yīng)用簡(jiǎn)介 6.5.3.1 偏微分方程求解程序概述,啟動(dòng)偏微分方程求解界面在 MATLAB 下鍵入 pdetool 該界面分為四個(gè)部分菜單系統(tǒng)工具欄集合編輯求解區(qū)域,6.5.3.2 偏微分方程求解區(qū)域繪

69、制,1）用工具欄中的橢圓、矩形等繪制一些區(qū)域。2）在集合編輯欄中修改其內(nèi)容。 如（R1＋E1＋E2）－E33）單擊工具欄中 按紐可得求解邊界。4）選擇Boundary-Remove All Subdomain Borders菜單項(xiàng)，消除相鄰區(qū)域中間的分隔線。5）單擊 按紐可將求解區(qū)域用三角形劃分成網(wǎng)格。可用 按紐加密。,,6.5.3.3 偏微分方程邊界條

70、件描述,選擇Boundary-Specify Boundary Conditions菜單狄利克雷條件，諾伊曼條件。,6.5.3.4 偏微分方程求解舉例,例：求解：1）繪制求解區(qū)域。2）描述邊界條件（Boundary-Specify Boundary Conditions）。3）選擇偏微分方程的類型。單擊工具欄中的PDE圖標(biāo)，在打開的新窗口選擇Hyperbolic選項(xiàng)，輸入?yún)?shù)c,a,f,d.4）求解。單擊工具欄中的等

71、號(hào)按鈕。,顯示：1）圖形顏色表示t=0時(shí)u(x,y)的函數(shù)值。2）單擊工具欄中的三維圖標(biāo)將打開一新的對(duì)話框，若再選擇Contour可繪制等值線，若選擇Arrows選項(xiàng)可繪制引力線。若單獨(dú)選擇Height（3d-plot)，則在另一窗口繪制出三維圖形。3）可在單擊三維圖標(biāo)打開的新對(duì)話框中，對(duì)Property欄目的各個(gè)項(xiàng)目重新選擇。4）可修改微分方程的邊界條件，重新求解。動(dòng)畫：1）Solve-Parameters對(duì)話框時(shí)間向量改

72、為0:0.1:2。2）三維圖標(biāo)打開的對(duì)話框中選擇Animation選項(xiàng)，單擊Options按紐設(shè)置播放速度。Plot-Export Movie 菜單。,6.5.3.5 函數(shù)參數(shù)的偏微分方程求解,例：(橢圓型）,求解：1）求解區(qū)域不變。2）描述邊界條件，u=0。3）選擇偏微分方程的類型。單擊工具欄中的PDE圖標(biāo)，在打開的新窗口選擇Elliptic選項(xiàng)，輸入?yún)?shù)c=1./sqrt(1+ux.^2+uy.^2), a=x.^2+y.^

73、2 , f=exp(-x.^2-y.^2).4)再打開Solve-Parameters對(duì)話框，選定Use nonlinear solve屬性（該屬性只適于橢圓性偏微分方程）5）求解。單擊工具欄中的等號(hào)按鈕。,例：橢圓微分方程的拉普拉斯形式：其自變量取值范圍D為： 邊界條件為：,1、調(diào)整坐標(biāo)范圍:options->Axes Limits2、設(shè)定矩形區(qū)域，雙擊確定調(diào)整。3、設(shè)置邊界條件：?jiǎn)螕暨吔绨粹o