偏微分方程數(shù)值解法地MATLAB源碼

1、實(shí)用標(biāo)準(zhǔn)原創(chuàng)偏微分方程數(shù)值解法的 MATLAEW碼【更新完畢】說(shuō)明：由于偏微分的程序都比較長(zhǎng)，比其他的算法稍復(fù)雜一些，所以另開(kāi)一貼，專門(mén)上傳偏微分的程序謝謝大家的支持！其他的數(shù)值算法見(jiàn)：.Announce/Announce.asp?BoardID=209&id=82450041、古典顯式格式求解拋物型偏微分方程(一維熱傳導(dǎo)方程)function U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%古典顯式格式求解拋物型偏微分方程%U x t=PDEParabolicClassicalExplicit(uX,uT,

2、phi,psi1,psi2,M,N,C)%程：u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT%初值條件：u(x,0)=phi(x)%&值條件：u(0,t)=psi1(t), u(uX,t)=psi2(t)%喻出參數(shù)：U -解矩陣，第一行表示初值，第一列和最后一列表示邊值，第二行表示第2層%x -空間變量%t -時(shí)間變量力輸入?yún)?shù)：uX -空間變量x的取值上限%uT -時(shí)間變量t的取值上限%phi -初值條件，定義為內(nèi)聯(lián)函數(shù)%psil -邊值條件，定義為內(nèi)聯(lián)函數(shù)%psi2 -邊值條件，定義為內(nèi)聯(lián)函數(shù)%M-沿x軸的等分區(qū)間數(shù)%N-沿t軸的等

3、分區(qū)間數(shù)%C-系數(shù)，默認(rèn)情況下 C=1%應(yīng)用舉例：%uX=1;uT=0.2;M=15;N=100;C=1;%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');%U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C);力役置參數(shù)C的默認(rèn)值if nargin=7C=1;end諭算步長(zhǎng)dx=uX/M;%x 的步長(zhǎng)dt=uT/N;%t 的步長(zhǎng)x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;% 步長(zhǎng)比

4、 r1=1-2*r;if r > 0.5disp('r > 0.5, 不穩(wěn)定') end%+算初值和邊值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i);endfor j=1:N+1U(1,j)=psi1(t(j);U(M+1,j)=psi2(t(j);end嗨層求解forj=1:Nfor i=2:MU(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j); endendU=U'%乍出圖形mesh(x,t,U);title('古典顯式格式，一維熱傳導(dǎo)方程的解的圖像')xlabel(

5、' 空間變量x')ylabel(' 時(shí)間變量t') zlabel('一維熱傳導(dǎo)方程的解U')return;古典顯式格式不穩(wěn)定情況文案大全-arsslklr隼一段 I古星rttlM- 一窕里恒寫(xiě)艮嗯kM的圖MJIM，一古典顯式格式穩(wěn)定情況Q DCLhdriUz cIlJL-LD二L iEh Te 琳 f Im«h iDo-k 3Hb干 1mMfr 14?西前詼籌冠苜石謨 1 tf同aIT ""* '空月磨置nli：3情戔*fT #-，g2、古典隱式格式求解拋物型偏微分方程(一維熱傳導(dǎo)方程)function U

6、x t=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C)“古典隱式格式求解拋物型偏微分方程%U x t=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C)%程：u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT“初值條件：u(x,0)=phi(x)“&值條件：u(0,t)=psi1(t), u(uX,t)=psi2(t)%喻出參數(shù)：U -解矩陣，第一行表示初值，第一列和最后一列表示邊值，第二行表示第2層%喻入?yún)?shù)x

7、-空間變量t -時(shí)間變量uX -空間變量x的取值上限%uT -時(shí)間變量t的取值上限%phi -初值條件，定義為內(nèi)聯(lián)函數(shù)%psi1 -邊值條件，定義為內(nèi)聯(lián)函數(shù)%psi2 -邊值條件，定義為內(nèi)聯(lián)函數(shù)%M -沿x軸的等分區(qū)間數(shù)%N -沿t軸的等分區(qū)間數(shù)%C -系數(shù)，默認(rèn)情況下 C=1%應(yīng)用舉例：%uX=1;uT=0.2;M=50;N=50;C=1;%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');%U x t=PDEParabolicClassicalImplicit(uX,uT,ph

8、i,psi1,psi2,M,N,C);力役置參數(shù)C的默認(rèn)值if nargin=7C=1;end%+算步長(zhǎng)dx=uX/M;%x 的步長(zhǎng) dt=uT/N;%t 的步長(zhǎng)x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;% 步長(zhǎng)比Diag=zeros(1,M-1);% 矩陣的對(duì)角線元素Low=zeros(1,M-2);%矩陣的下對(duì)角線元素 Up=zeros(1,M-2);%矩陣的上對(duì)角線元素 for i=1:M-2Diag(i)=1+2*r;Low(i)=-r;Up(i)=-r;endDiag(M-1)=1+2*r;%+算初值和邊值U=zeros(M+1,N+1);for i=1:

9、M+1U(i,1)=phi(x(i);endfor j=1:N+1U(1,j)=psi1(t(j);U(M+1,j戶psi2(t(j);end%逐層求解，需要使用追趕法(調(diào)用函數(shù)EqtsForwardAndBackward )forj=1:Nb1=zeros(M-1,1);b1(1)=r*U(1,j+1);b1(M-1)=r*U(M+1,j+1);b=U(2:M,j)+b1;U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);endU=U'%乍出圖形mesh(x,t,U);title('古典隱式格式，一維熱傳導(dǎo)方程的解的圖像'

10、;)xlabel(' 空間變量x')ylabel(' 時(shí)間變量t')zlabel(' 一維熱傳導(dǎo)方程的解U')return;此算法需要使用追趕法求解三對(duì)角線性方程組，這個(gè)算法在上一篇帖子中已經(jīng)給出，為了方便，再給出來(lái)追趕法解三對(duì)角線性方程組function x=EqtsForwardAndBackward(L,D,U,b)%追趕法求解三對(duì)角線性方程組Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三對(duì)角線性方程組的解%L:三對(duì)角矩陣的下對(duì)角線，行向量%口三對(duì)角矩陣的對(duì)角線，行向量%U三對(duì)角矩陣的上對(duì)角線，行向量%

11、b:線性方程組Ax=b中的b,列向量%應(yīng)用舉例：%L=-1 -2 -3;D=2 3 4 5;U=-1 -2 -3;b=6 1 -2 1'%x=EqtsForwardAndBackward(L,D,U,b)%檢查參數(shù)的輸入是否正確n=length(D);m=length(b);n1=length(L);n2=length(U);if n-n1 = 1 | n-n2 = 1 11n = m disp('輸入?yún)?shù)有誤！) x=''return;end%追的過(guò)程for i=2:nL(i-1)=L(i-1)/D(i-1);D(i)=D(i)-L(i-1)*U(i-1);e

12、ndx=zeros(n,1);x(1)=b(1);for i=2:nx(i)=b(i)-L(i-1)*x(i-1);end%趕的過(guò)程x(n)=x(n)/D(n);for i=n-1:-1:1x(i)=(x(i)-U(i)*x(i+1)/D(i);endreturn;古典隱式格式在以后的程序中，我們都取 C=1,不再作為一個(gè)輸入?yún)?shù)處理3、Crank-Nicolson隱式格式求解拋物型偏微分方程需要調(diào)用追趕法的程序function U x t=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)%Crank-Nicolson隱式格式求解拋物型偏微分方程%U x t=PD

13、EParabolicCN(uX,uT,phi,psi1,psi2,M,N) %程：u_t=u_xx 0 <= x <= uX,0 <= t <= uT%初值條件：u(x,0)=phi(x)%i值條件：u(0,t)=psi1(t), u(uX,t)=psi2(t)%喻出參數(shù)：U -解矩陣，第一行表示初值，第一列和最后一列表示邊值，第二行表示第2層%x -空間變量%t -時(shí)間變量力輸入?yún)?shù)：uX -空間變量x的取值上限%uT -時(shí)間變量t的取值上限%phi -初值條件，定義為內(nèi)聯(lián)函數(shù)%psil -邊值條件，定義為內(nèi)聯(lián)函數(shù)%psi2 -邊值條件，定義為內(nèi)聯(lián)函數(shù)%M -沿x軸的

14、等分區(qū)間數(shù)%N -沿t軸的等分區(qū)間數(shù)%應(yīng)用舉例：%uX=1;uT=0.2;M=50;N=50;%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');%U x t=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N);%+算步長(zhǎng)dx=uX/M;%x 的步長(zhǎng)dt=uT/N;%t 的步長(zhǎng)x=(0:M)*dx;t=(0:N)*dt;r=dt/dx/dx;% 步長(zhǎng)比Diag=zeros(1,M-1);% 矩陣的對(duì)角線元素Low=zeros(1,M-2);%矩陣的下對(duì)角

15、線元素Up=zeros(1,M-2);%矩陣的上對(duì)角線元素for i=1:M-2Diag(i)=1+r;Low(i)=-r/2;Up(i)=-r/2;endDiag(M-1)=1+r;%+算初值和邊值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i);endfor j=1:N+1U(1,j)=psi1(t(j);U(M+1,j)=psi2(t(j);endB=zeros(M-1,M-1);for i=1:M-2B(i,i)=1-r;B(i,i+1)=r/2;B(i+1,i)=r/2;endB(M-1,M-1)=1-r;%逐層求解，需要使用追趕法(調(diào)用函數(shù) E

16、qtsForwardAndBackward )forj=1:Nb1=zeros(M-1,1);b1(1)=r*(U(1,j+1)+U(1,j)/2;b1(M-1)=r*(U(M+1,j+1)+U(M+1,j)/2;b=B*U(2:M,j)+b1;U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);endU=U'%乍出圖形mesh(x,t,U);title('Crank-Nicolson隱式格式，一維熱傳導(dǎo)方程的解的圖像')xlabel(' 空間變量x')ylabel(' 時(shí)間變量t')zlab

17、el(' 一維熱傳導(dǎo)方程的解U')return;隱式格式Crank-Nicolson£»* H mis m懇立電It 一 一事羈悔K弓g矍w耀.4、正方形區(qū)域 Laplace方程Diriclet 問(wèn)題的求解需要調(diào)用Jacobi迭代法和Guass-Seidel迭代法求解線性方程組function U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%SE方形區(qū)域Laplace方程的Diriclet邊值問(wèn)題的差分求解%t程序需要調(diào)用Jacobi迭代法或者Guass-Seidel

18、迭代法求解線性方程組%U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%5"程：u_xx+u_yy=0 0<=x,y<=ub%& 值條件：u(0,y)=phi1(y)%u(ub,y)=phi2(y)%u(x,0)=psi1(x)%u(x,ub)=psi2(x)%喻出參數(shù)：U -解矩陣，第一行表示 y=0時(shí)的值，第二行表示第 y=h時(shí)的值%x -橫坐標(biāo)%y -縱坐標(biāo)%輸入?yún)?shù)：ub -變量邊界值的上限%phi1,phi2,psi1,psi2 -邊界函數(shù)，定義為內(nèi)聯(lián)函數(shù)%M -橫

19、縱坐標(biāo)的等分區(qū)間數(shù)%type -求解差分方程的迭代格式，若 type='Jacobi',采用Jacobi迭代格式%若type='GS',采用Guass-Seidel迭代格式。默認(rèn)情況下，type='GS'%應(yīng)用舉例：%ub=4;M=20;%phi1=inline('y*(4-y)');phi2=inline('0');psi1=inline('sin(pi*x/4)');psi2=inline('0');%U x y=PDEEllipseSquareLaplaceDirichlet(

20、ub,phi1,phi2,psi1,psi2,M,'GS');if nargin=6type='GS'end唯長(zhǎng)h=ub/M;喊縱坐標(biāo)x=(0:M)*h;y=(0:M)*h;播分格式的矩陣形式 AU=K%勾造矢!陣AM2=(M-1)A2;A=zeros(M2);for i=1:M2A(i,i)=4;endfor i=1:M2-1if mod(i,M-1)=0A(i,i+1)=-1;A(i+1,i)=-1;endendfor i=1:M2-M+1A(i,i+M-1)=-1;A(i+M-1,i)=-1;endU=zeros(M+1);%i值條件for i=1:M+1

21、U(i,1)=psi1(i-1)*h);U(i,M+1)=psi2(i-1)*h);U(1,i)=phi1(i-1)*h);U(M+1,i)=phi2(i-1)*h);end%勾造KK=zeros(M2,1);for i=1:M-1K(i)=U(i+1,1);K(M2-i+1)=U(i+1,M+1);endK(1)=K(1)+U(1,2);K(M-1)=K(M-1)+U(M+1,2);K(M2-M+2)=K(M2-M+2)+U(1,M);K(M2)=K(M2)+U(M+1,M);for i=2:M-2K(M-1)*(i-1)+1)=U(1,i+1);K(M-1)*i)=U(M+1,i+1);e

22、ndx0=ones(M2,1);switch type%調(diào)用Guass-Seidel迭代法求解線性方程組AU=Kcase 'Jacobi'X=EqtsJacobi(A,K,x0);%調(diào)用Guass-Seidel迭代法求解線性方程組AU=Kcase 'GS'X=EqtsGS(A,K,x0);otherwisedisp('差分格式類型輸入錯(cuò)誤')return;end%巴求解結(jié)果化成矩陣型式for i=2:Mfor j=2:MU(j,i)=X(j-1+(M-1)*(i-2);endendU=U'%乍出圖形mesh(x,y,U);title(&#

23、39; 五點(diǎn)差分格式Laplace方程Diriclet問(wèn)題的解的圖像')xlabel('x')ylabel('y')zlabel('Laplace 方程 Diriclet 問(wèn)題的解 U')return;河正正方形區(qū)域JIMLaplace方程五點(diǎn)差分格式i£4e 出f levin 匕卜 EMilAfip 貼T»r %.酋荀箭甲TRW孱罰a a-E n14 - -5/ l：L. *而 H國(guó)*5、一階雙曲型方程的差分方法function U x t=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi

24、2,type)%廠階雙曲型方程的差分格式%U x t=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)%f 程:u_t+C*u_x=00 <= t <= uT, 0 <= x <= uX%初值條件：u(x,0)=phi(x)%喻出參數(shù):U -解矩陣，第一行表示初值，第二行表示第2個(gè)時(shí)間層%x -橫坐標(biāo)%t -縱坐標(biāo)，時(shí)間%輸入?yún)?shù)：uX -變量x的上界%uT -變量t的上界%M -變量x的等分區(qū)間數(shù)%N -變量t的等分區(qū)間數(shù)%C -系數(shù)%phi-初值條件函數(shù)，定義為內(nèi)聯(lián)函數(shù)%psi1,psi2 - 邊值條件函數(shù)，定義為內(nèi)聯(lián)函數(shù)%

25、type -差分格式，從下列值中選取%-type='LaxFriedrichs' ,采用 Lax-Friedrichs 差分格式求解%-type='CourantIsaacsonRees' ,采用 Courant-Isaacson-Rees 差分格式求解%-type='LeapFrog',采用Leap-Frog (蛙跳)差分格式求解%-type='LaxWendroff',采用 Lax-Wendroff 差分格式求解%-type='CrankNicolson' ,采用 Crank-Nicolson 差分格式求解，此

26、格式需調(diào)用追趕法%求解三對(duì)角線性方程組%h=uX/M;%變量x的步長(zhǎng)k=uT/N;%變量t的步長(zhǎng) r=k/h;%步長(zhǎng)比x=(0:M)*h;t=(0:N)*k;U=zeros(M+1,N+1);力初值條件for i=1:M+1U(i,1)=phi(x(i);end%邊值條件for j=1:N+1U(1,j)=psi1(t(j);U(M+1,j)=psi2(t(j);%U(1,j)=NaN;%U(M+1,j)=NaN;endswitch type %Lax-Friedrichs 差分格式case 'LaxFriedrichs'if abs(C*r)>1disp('|C

27、*r|>1, Lax-Friedrichs差分格式不穩(wěn)定！')end%逐層求解forj=1:Nfor i=2:MU(i,j+1)=(U(i+1,j)+U(i-1,j)/2-C*r*(U(i+1,j)-U(i-1,j)/2;endend%Courant-Isaacson-Rees 差分格式 case 'CourantlsaacsonRees' if C<0 disp('C<0 ,采用前差公式') if C*r<-1disp('Courant-Isaacson-Lees差分格式不穩(wěn)定！')end %逐層求解 forj=

28、1:N for i=2:MU(i,j+1)=(1+C*r)*U(i,j)-C*r*U(i+1,j);endendelsedisp('C>0 ,采用后差公式') if C*r>1disp('Courant-Isaacson-Lees差分格式不穩(wěn)定！')end %逐層求解 forj=1:N for i=2:MU(i,j+1)=C*r*U(i-1,j)+(1-C*r)*U(i,j);endendend%Leap-Frog (蛙跳)差分格式 case 'LeapFrog'phi2=input('請(qǐng)輸入第二層初值條件函數(shù)：psi2=');if abs(C*r)>1disp('|C*r|>1, Leap-Frog 差分格式不穩(wěn)定！')end %第二層初值條件 for i=1:M+1U(i,2)=phi2(x(i);end%逐層求解for j=2:Nfor i=2:MU(i,j+1