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

1、.原創(chuàng)偏微分方程數(shù)值解法的MATLAB源碼【更新完畢】說明：由于偏微分的程序都比較長，比其他的算法稍復(fù)雜一些，所以另開一貼，專門上傳偏微分的程序謝謝大家的支持！其他的數(shù)值算法見：./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,phi

2、,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 -時間變量%輸入?yún)?shù)：uX -空間變量x的取值上限%       &

3、#160; uT -時間變量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ù)%       

4、  C -系數(shù)，默認情況下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è)置參數(shù)C的默認值if nargin=7    C=1;end%計算步長dx=uX/M;%x的步長dt=uT/N;%t的步長x=(0:M)*dx;t=(0:N)*dt;r=C

5、*dt/dx/dx;%步長比r1=1-2*r;if r > 0.5    disp('r > 0.5,不穩(wěn)定')end%計算初值和邊值U=zeros(M+1,N+1);for i=1:M+1    U(i,1)=phi(x(i);endfor j=1:N+1    U(1,j)=psi1(t(j);    U(M+1,j)=psi2(t(j);end%逐層求解for j=1:N    for i=2:M        

6、U(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('空間變量 x')ylabel('時間變量 t')zlabel('一維熱傳導(dǎo)方程的解 U')return;古典顯式格式不穩(wěn)定情況古典顯式格式穩(wěn)定情況2、古典隱式格式求解拋物型偏微分方程（一維熱傳導(dǎo)方程）function U x t=PDEParabolicClassicalImplicit(u

7、X,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層%         x -空間變量% &

8、#160;       t -時間變量%輸入?yún)?shù)：uX -空間變量x的取值上限%         uT -時間變量t的取值上限%         phi -初值條件，定義為內(nèi)聯(lián)函數(shù)%         psi1 -邊值條件，定義為內(nèi)聯(lián)函數(shù)%         psi2 -邊值條件，定義為內(nèi)聯(lián)函數(shù)%        

9、 M -沿x軸的等分區(qū)間數(shù)%         N -沿t軸的等分區(qū)間數(shù)%         C -系數(shù)，默認情況下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,phi,psi1,psi2,M,N,C);%設(shè)置參數(shù)C的默認值i

10、f nargin=7    C=1;end%計算步長dx=uX/M;%x的步長dt=uT/N;%t的步長x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;%步長比Diag=zeros(1,M-1);%矩陣的對角線元素Low=zeros(1,M-2);%矩陣的下對角線元素Up=zeros(1,M-2);%矩陣的上對角線元素for i=1:M-2    Diag(i)=1+2*r;    Low(i)=-r;    Up(i)=-r;endDiag(M-1)=1+2*r;%計算初值和邊值U=zero

11、s(M+1,N+1);for i=1:M+1    U(i,1)=phi(x(i);endfor j=1:N+1    U(1,j)=psi1(t(j);    U(M+1,j)=psi2(t(j);end%逐層求解，需要使用追趕法（調(diào)用函數(shù)EqtsForwardAndBackward）for j=1:N    b1=zeros(M-1,1);    b1(1)=r*U(1,j+1);    b1(M-1)=r*U(M+1,j+1);    b=U(2:M

12、,j)+b1;    U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);endU=U'%作出圖形mesh(x,t,U);title('古典隱式格式，一維熱傳導(dǎo)方程的解的圖像')xlabel('空間變量 x')ylabel('時間變量 t')zlabel('一維熱傳導(dǎo)方程的解 U')return;此算法需要使用追趕法求解三對角線性方程組，這個算法在上一篇帖子中已經(jīng)給出，為了方便，再給出來追趕法解三對角線性方程組function x=EqtsForwardAn

13、dBackward(L,D,U,b)%追趕法求解三對角線性方程組Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三對角線性方程組的解%L:三對角矩陣的下對角線，行向量%D:三對角矩陣的對角線，行向量%U:三對角矩陣的上對角線，行向量%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

14、);if n-n1 = 1 | n-n2 = 1 | n = m    disp('輸入?yún)?shù)有誤！')    x=' '    return;end%追的過程for i=2:n    L(i-1)=L(i-1)/D(i-1);    D(i)=D(i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);for i=2:n    x(i)=b(i)-L(i-1)*x(i-1);end%趕的過程x(n)=x(n)/D(

15、n);for i=n-1:-1:1    x(i)=(x(i)-U(i)*x(i+1)/D(i);endreturn;古典隱式格式在以后的程序中，我們都取C=1，不再作為一個輸入?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=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)%方程：u_t=u_xx 0 <= x <=

16、 uX,0 <= t <= uT%初值條件：u(x,0)=phi(x)%邊值條件：u(0,t)=psi1(t), u(uX,t)=psi2(t)%輸出參數(shù)：U -解矩陣，第一行表示初值，第一列和最后一列表示邊值，第二行表示第2層%         x -空間變量%         t -時間變量%輸入?yún)?shù)：uX -空間變量x的取值上限%         uT -時間變量t的取值上限%       

17、;  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ù)%應(yīng)用舉例：%uX=1;uT=0.2;M=50;N=50;%phi=inline('sin(pi*x)');psi1=inline(&

18、#39;0');psi2=inline('0');%U x t=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N);%計算步長dx=uX/M;%x的步長dt=uT/N;%t的步長x=(0:M)*dx;t=(0:N)*dt;r=dt/dx/dx;%步長比Diag=zeros(1,M-1);%矩陣的對角線元素Low=zeros(1,M-2);%矩陣的下對角線元素Up=zeros(1,M-2);%矩陣的上對角線元素for i=1:M-2    Diag(i)=1+r;    Low(i)=-r/2; 

19、;   Up(i)=-r/2;endDiag(M-1)=1+r;%計算初值和邊值U=zeros(M+1,N+1);for i=1:M+1    U(i,1)=phi(x(i);endfor j=1:N+1    U(1,j)=psi1(t(j);    U(M+1,j)=psi2(t(j);endB=zeros(M-1,M-1);for i=1:M-2    B(i,i)=1-r;    B(i,i+1)=r/2;    B(i+1,i)=r/2;endB(M-1,

20、M-1)=1-r;%逐層求解，需要使用追趕法（調(diào)用函數(shù)EqtsForwardAndBackward）for j=1:N    b1=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);ti

21、tle('Crank-Nicolson隱式格式，一維熱傳導(dǎo)方程的解的圖像')xlabel('空間變量 x')ylabel('時間變量 t')zlabel('一維熱傳導(dǎo)方程的解 U')return;Crank-Nicolson隱式格式4、正方形區(qū)域Laplace方程Diriclet問題的求解需要調(diào)用Jacobi迭代法和Guass-Seidel迭代法求解線性方程組function U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%正方形區(qū)域Lap

22、lace方程的Diriclet邊值問題的差分求解%此程序需要調(diào)用Jacobi迭代法或者Guass-Seidel迭代法求解線性方程組%U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%方程：u_xx+u_yy=0  0<=x,y<=ub%邊值條件：u(0,y)=phi1(y)%         u(ub,y)=phi2(y)%         u(x,0)=psi1(x)%  

23、       u(x,ub)=psi2(x)%輸出參數(shù)：U -解矩陣，第一行表示y=0時的值，第二行表示第y=h時的值%         x -橫坐標%         y -縱坐標%輸入?yún)?shù)：ub -變量邊界值的上限%         phi1,phi2,psi1,psi2 -邊界函數(shù)，定義為內(nèi)聯(lián)函數(shù)%         M -橫縱坐標的等分區(qū)間數(shù)%

24、60;       type -求解差分方程的迭代格式，若type='Jacobi'，采用Jacobi迭代格式%               若type='GS'，采用Guass-Seidel迭代格式。默認情況下，type='GS'%應(yīng)用舉例：%ub=4;M=20;%phi1=inline('y*(4-y)');phi2=inline('0');psi1=inline('sin(

25、pi*x/4)');psi2=inline('0');%U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,'GS');if nargin=6    type='GS'end%步長h=ub/M;%橫縱坐標x=(0:M)*h;y=(0:M)*h;%差分格式的矩陣形式AU=K%構(gòu)造矩陣AM2=(M-1)2;A=zeros(M2);for i=1:M2    A(i,i)=4;endfor i=1:M2-1   

26、; if mod(i,M-1)=0        A(i,i+1)=-1;        A(i+1,i)=-1;    endendfor i=1:M2-M+1    A(i,i+M-1)=-1;    A(i+M-1,i)=-1;endU=zeros(M+1);%邊值條件for i=1:M+1    U(i,1)=psi1(i-1)*h);    U(i,M+1)=psi2(i-1)*h);&

27、#160;   U(1,i)=phi1(i-1)*h);    U(M+1,i)=phi2(i-1)*h);end%構(gòu)造KK=zeros(M2,1);for i=1:M-1    K(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-2    K(M-1)*(i-1)+1)=U(1,

28、i+1);    K(M-1)*i)=U(M+1,i+1);endx0=ones(M2,1);switch type    %調(diào)用Guass-Seidel迭代法求解線性方程組AU=K    case 'Jacobi'        X=EqtsJacobi(A,K,x0);    %調(diào)用Guass-Seidel迭代法求解線性方程組AU=K    case 'GS'      

29、  X=EqtsGS(A,K,x0);    otherwise        disp('差分格式類型輸入錯誤')        return;end%把求解結(jié)果化成矩陣型式for i=2:M    for j=2:M        U(j,i)=X(j-1+(M-1)*(i-2);    endendU=U'%作出圖形mesh(x,y,U);titl

30、e('五點差分格式Laplace方程Diriclet問題的解的圖像')xlabel('x')ylabel('y')zlabel('Laplace方程Diriclet問題的解 U')return;正方形區(qū)域Laplace方程五點差分格式5、一階雙曲型方程的差分方法function U x t=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)%一階雙曲型方程的差分格式%U x t=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)%方程:u_t+C*u_

31、x=0  0 <= t <= uT, 0 <= x <= uX%初值條件:u(x,0)=phi(x)%輸出參數(shù):U -解矩陣，第一行表示初值，第二行表示第2個時間層%        x -橫坐標%        t -縱坐標，時間%輸入?yún)?shù):uX -變量x的上界%        uT -變量t的上界%        M -變量x的等分區(qū)間數(shù)%        N -

32、變量t的等分區(qū)間數(shù)%        C -系數(shù)%        phi -初值條件函數(shù)，定義為內(nèi)聯(lián)函數(shù)%        psi1,psi2 -邊值條件函數(shù)，定義為內(nèi)聯(lián)函數(shù)%        type -差分格式，從下列值中選取%             -type='LaxFriedrichs'，采用Lax-Friedrichs差分格式求解%

33、            -type='CourantIsaacsonRees'，采用Courant-Isaacson-Rees差分格式求解%             -type='LeapFrog'，采用Leap-Frog（蛙跳）差分格式求解%             -type='LaxWendroff'，采

34、用Lax-Wendroff差分格式求解%             -type='CrankNicolson'，采用Crank-Nicolson差分格式求解，此格式需調(diào)用追趕法%              求解三對角線性方程組%h=uX/M;%變量x的步長k=uT/N;%變量t的步長r=k/h;%步長比x=(0:M)*h;t=(0:N)*k;U=zeros(M+1,N+1);%初值條件for i=1:M+1  &

35、#160; U(i,1)=phi(x(i);end%邊值條件for j=1:N+1    U(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)>1

36、60;           disp('|C*r|>1，Lax-Friedrichs差分格式不穩(wěn)定！')        end        %逐層求解        for j=1:N            for i=2:M        &#

37、160;       U(i,j+1)=(U(i+1,j)+U(i-1,j)/2-C*r*(U(i+1,j)-U(i-1,j)/2;            end        end        %Courant-Isaacson-Rees差分格式    case 'CourantIsaacsonRees'      

38、  if C<0            disp('C<0，采用前差公式')            if C*r<-1                disp('Courant-Isaacson-Lees差分格式不穩(wěn)定！')         

39、;   end            %逐層求解            for j=1:N                for i=2:M                    U(i,j+1)=(1+C*r)*U(i

40、,j)-C*r*U(i+1,j);                end            end        else            disp('C>0，采用后差公式')            if C

41、*r>1                disp('Courant-Isaacson-Lees差分格式不穩(wěn)定！')            end            %逐層求解            for j=1:N     

42、0;          for i=2:M                    U(i,j+1)=C*r*U(i-1,j)+(1-C*r)*U(i,j);                end            end  &

43、#160;     end            %Leap-Frog（蛙跳）差分格式    case 'LeapFrog'        phi2=input('請輸入第二層初值條件函數(shù)：psi2=');        if abs(C*r)>1           

44、disp('|C*r|>1，Leap-Frog差分格式不穩(wěn)定！')        end        %第二層初值條件        for i=1:M+1            U(i,2)=phi2(x(i);        end        %逐

45、層求解        for j=2:N            for i=2:M                U(i,j+1)=U(i,j-1)-C*r*(U(i+1,j)-U(i-1,j);            end        end 

46、           %Lax-Wendroff差分格式    case 'LaxWendroff'        if abs(C*r)>1            disp('|C*r|>1，Lax-Wendroff差分格式不穩(wěn)定！')        end     

47、60;  %逐層求解        for j=1:N            for i=2:M                U(i,j+1)=U(i,j)-C*r*(U(i+1,j)-U(i-1,j)/2+C2*r2*(U(i+1,j)-2*U(i,j)+U(i-1,j)/2;            end        end            %Crank-Nicolson隱式差分格式，需調(diào)用追趕法求解三對角線性方程組的算法    case 'CrankNicolson'        Diag=zeros(1,M-1);%矩陣的對角線元素     &