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

《偏微分方程數(shù)值解法的MATLAB源碼》由會(huì)員分享，可在線閱讀，更多相關(guān)《偏微分方程數(shù)值解法的MATLAB源碼（16頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

1、精選優(yōu)質(zhì)文檔-----傾情為你奉上 [ 原創(chuàng)]偏微分方程數(shù)值解法的MATLAB源碼【更新完畢】 說(shuō)明：由于偏微分的程序都比較長(zhǎng)，比其他的算法稍復(fù)雜一些，所以另開(kāi)一貼，專門上傳偏微分的程序 謝謝大家的支持！ 其他的數(shù)值算法見(jiàn)： 1、古典顯式格式求解拋物型偏微分方程（一維熱傳導(dǎo)方程） function [U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典顯式格式求解拋物型偏微分方程 %[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,

2、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ù) % ? ?? ?? psi1 -邊值條件，定義為內(nèi)聯(lián)函數(shù) % ? ??

3、 ?? 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=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的默認(rèn)值 if nargin==7 ? ? C=1; end %計(jì)

4、算步長(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)比 r1=1-2*r; if r > 0.5 ? ? disp('r > 0.5,不穩(wěn)定') end %計(jì)算初值和邊值 U=zeros(M+1,N+1); for i=1:M+1 ? ? 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)); end %逐層求解 for j=1:N ? ? for i=2:M ? ? ?? ?U(i,

5、j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j); ? ? end end U=U'; %作出圖形 mesh(x,t,U); title('古典顯式格式，一維熱傳導(dǎo)方程的解的圖像') xlabel('空間變量 x') ylabel('時(shí)間變量 t') zlabel('一維熱傳導(dǎo)方程的解 U') return; 古典顯式格式不穩(wěn)定情況 古典顯式格式穩(wěn)定情況 2、古典隱式格式求解拋物型偏微分方程（一維熱傳導(dǎo)方程） function [U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,

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

7、 ?? 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]=PDEParabolicClassi

8、calImplicit(uX,uT,phi,psi1,psi2,M,N,C); %設(shè)置參數(shù)C的默認(rèn)值 if nargin==7 ? ? C=1; end %計(jì)算步長(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-2 ? ? Diag(i)=1+2*r; ? ? Low(i)=-r; ? ? Up(i)=-r

9、; end Diag(M-1)=1+2*r; %計(jì)算初值和邊值 U=zeros(M+1,N+1); for i=1:M+1 ? ? 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)); 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,j)+b1; ? ? U(

10、2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b); end U=U'; %作出圖形 mesh(x,t,U); title('古典隱式格式，一維熱傳導(dǎo)方程的解的圖像') 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=

11、b %x=EqtsForwardAndBackward(L,D,U,b) %x:三對(duì)角線性方程組的解 %L:三對(duì)角矩陣的下對(duì)角線，行向量 %D:三對(duì)角矩陣的對(duì)角線，行向量 %U:三對(duì)角矩陣的上對(duì)角線，行向量 %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

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

13、 古典隱式格式 在以后的程序中，我們都取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]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N) % %方程：u_t=u_xx 0 <= x <= uX,0 <= t <= uT %初值條件：u(x,0)=phi(x) %邊值條件：u(0,t)=psi1

14、(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ù) % ? ?? ?? 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=

15、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); %計(jì)算步長(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ì)角線元素 Up=zeros(1,M-2);%矩陣的上對(duì)角線元素 for i=1:M-2 ?

16、 ? Diag(i)=1+r; ? ? Low(i)=-r/2; ? ? Up(i)=-r/2; end Diag(M-1)=1+r; %計(jì)算初值和邊值 U=zeros(M+1,N+1); for i=1:M+1 ? ? 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)); end B=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; end B(M-1,M-1)=1-r;

17、 %逐層求解，需要使用追趕法（調(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); end U=U'; %作出圖形 mesh(x,t,U); title('Crank-Nicolson隱式格式，一維熱傳導(dǎo)方程的解的圖像') xlabel

18、('空間變量 x') ylabel('時(shí)間變量 t') zlabel('一維熱傳導(dǎo)方程的解 U') return; Crank-Nicolson隱式格式 4、正方形區(qū)域Laplace方程Diriclet問(wèn)題的求解 需要調(diào)用Jacobi迭代法和Guass-Seidel迭代法求解線性方程組 function [U x y]=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type) %正方形區(qū)域Laplace方程的Diriclet邊值問(wèn)題的差分求解 %此程序需要調(diào)用Jacobi迭代法或者Guass-Seidel

19、迭代法求解線性方程組 %[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) % ? ?? ?? u(x,ub)=psi2(x) % %輸出參數(shù)：U -解矩陣，第一行表示y=0時(shí)的值，第二行表示第y=h時(shí)的值…… % ? ?? ?? x -橫坐標(biāo) % ? ?? ?? y -縱坐標(biāo) %輸入?yún)?shù)：ub -變量

20、邊界值的上限 % ? ?? ?? phi1,phi2,psi1,psi2 -邊界函數(shù)，定義為內(nèi)聯(lián)函數(shù) % ? ?? ?? M -橫縱坐標(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');

21、%[U x y]=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,'GS'); if nargin==6 ? ? type='GS'; end %步長(zhǎng) h=ub/M; %橫縱坐標(biāo) x=(0:M)*h; y=(0:M)*h; %差分格式的矩陣形式AU=K %構(gòu)造矩陣A M2=(M-1)^2; A=zeros(M2); for i=1:M2 ? ? A(i,i)=4; end for i=1:M2-1 ? ? if mod(i,M-1)~=0 ? ? ?? ?A(i,i+1)=-1; ? ? ?? ?A(i+1,i)=-1

22、; ? ? end end for i=1:M2-M+1 ? ? A(i,i+M-1)=-1; ? ? A(i+M-1,i)=-1; end U=zeros(M+1); %邊值條件 for i=1:M+1 ? ? 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 %構(gòu)造K K=zeros(M2,1); for i=1:M-1 ? ? K(i)=U(i+1,1); ? ? K(M2-i+1)=U(i+1,M+1); end K(

23、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,i+1); ? ? K((M-1)*i)=U(M+1,i+1); end x0=ones(M2,1); switch type ? ? %調(diào)用Guass-Seidel迭代法求解線性方程組AU=K ? ? case 'Jacobi' ? ? ?? ?X=EqtsJacobi(A,K,x0); ? ? %調(diào)用Guass-Seidel迭代法求解線

24、性方程組AU=K ? ? case 'GS' ? ? ?? ?X=EqtsGS(A,K,x0); ? ? otherwise ? ? ?? ?disp('差分格式類型輸入錯(cuò)誤') ? ? ?? ?return; end %把求解結(jié)果化成矩陣型式 for i=2:M ? ? for j=2:M ? ? ?? ?U(j,i)=X(j-1+(M-1)*(i-2)); ? ? end end U=U'; %作出圖形 mesh(x,y,U); title('五點(diǎn)差分格式Laplace方程Diriclet問(wèn)題的解的圖像') xlabel('x') ylabel('y') zlabel('Laplace

25、方程Diriclet問(wèn)題的解 U') return; 正方形區(qū)域Laplace方程五點(diǎn)差分格式 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_x=0 ?0 <= t <= uT, 0 <= x <= uX %初值條件:u(x,0)=phi(x) % %輸出參數(shù):U -解矩陣，第一行表示初值，第二行表示第

26、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ù) % ? ?? ? type -差分格式，從下列值中選取 % ? ?? ?? ?? ?-type='LaxFriedrichs'，采用Lax-Friedrichs差分格式求解 % ? ?? ??

27、?? ?-type='CourantIsaacsonRees'，采用Courant-Isaacson-Rees差分格式求解 % ? ?? ?? ?? ?-type='LeapFrog'，采用Leap-Frog（蛙跳）差分格式求解 % ? ?? ?? ?? ?-type='LaxWendroff'，采用Lax-Wendroff差分格式求解 % ? ?? ?? ?? ?-type='CrankNicolson'，采用Crank-Nicolson差分格式求解，此格式需調(diào)用追趕法 % ? ?? ?? ?? ? 求解三對(duì)角線性方程組 % h=uX/M;%變量x的步長(zhǎng) k=uT/N;%變量t的步長(zhǎng) r=

28、k/h;%步長(zhǎng)比 x=(0:M)*h; t=(0:N)*k; U=zeros(M+1,N+1); %初值條件 for i=1:M+1 ? ? 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; end switch type ? ? %Lax-Friedrichs差分格式 ? ? case 'LaxFriedrichs' ? ? ?? ?if abs(C*r)>1 ? ? ?? ?? ? dis

29、p('|C*r|>1，Lax-Friedrichs差分格式不穩(wěn)定！') ? ? ?? ?end ? ? ?? ?%逐層求解 ? ? ?? ?for j=1:N ? ? ?? ?? ? for i=2:M ? ? ?? ?? ?? ?? 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' ? ? ?? ?if C<0 ? ? ?? ?? ? di

30、sp('C<0，采用前差公式') ? ? ?? ?? ? if C*r<-1 ? ? ?? ?? ?? ?? disp('Courant-Isaacson-Lees差分格式不穩(wěn)定！') ? ? ?? ?? ? end ? ? ?? ?? ? %逐層求解 ? ? ?? ?? ? for j=1:N ? ? ?? ?? ?? ?? for i=2:M ? ? ?? ?? ?? ?? ?? ?U(i,j+1)=(1+C*r)*U(i,j)-C*r*U(i+1,j); ? ? ?? ?? ?? ?? end ? ? ?? ?? ? end ? ? ?? ?else ? ? ?? ?? ? disp('

31、C>0，采用后差公式') ? ? ?? ?? ? if C*r>1 ? ? ?? ?? ?? ?? disp('Courant-Isaacson-Lees差分格式不穩(wěn)定！') ? ? ?? ?? ? end ? ? ?? ?? ? %逐層求解 ? ? ?? ?? ? for j=1:N ? ? ?? ?? ?? ?? for i=2:M ? ? ?? ?? ?? ?? ?? ?U(i,j+1)=C*r*U(i-1,j)+(1-C*r)*U(i,j); ? ? ?? ?? ?? ?? end ? ? ?? ?? ? end ? ? ?? ?end ? ? ?? ? ? ? %Leap-Frog（

32、蛙跳）差分格式 ? ? case 'LeapFrog' ? ? ?? ?phi2=input('請(qǐng)輸入第二層初值條件函數(shù)：psi2='); ? ? ?? ?if abs(C*r)>1 ? ? ?? ?? ? disp('|C*r|>1，Leap-Frog差分格式不穩(wěn)定！') ? ? ?? ?end ? ? ?? ?%第二層初值條件 ? ? ?? ?for i=1:M+1 ? ? ?? ?? ? U(i,2)=phi2(x(i)); ? ? ?? ?end ? ? ?? ?%逐層求解 ? ? ?? ?for j=2:N ? ? ?? ?? ? for i=2:M ? ? ?? ?? ?? ??

33、U(i,j+1)=U(i,j-1)-C*r*(U(i+1,j)-U(i-1,j)); ? ? ?? ?? ? end ? ? ?? ?end ? ? ?? ? ? ? %Lax-Wendroff差分格式 ? ? case 'LaxWendroff' ? ? ?? ?if abs(C*r)>1 ? ? ?? ?? ? disp('|C*r|>1，Lax-Wendroff差分格式不穩(wěn)定！') ? ? ?? ?end ? ? ?? ?%逐層求解 ? ? ?? ?for j=1:N ? ? ?? ?? ? for i=2:M ? ? ?? ?? ?? ?? U(i,j+1)=U(i,j)-C*r*(U

34、(i+1,j)-U(i-1,j))/2+C^2*r^2*(U(i+1,j)-2*U(i,j)+U(i-1,j))/2; ? ? ?? ?? ? end ? ? ?? ?end ? ? ?? ? ? ? %Crank-Nicolson隱式差分格式，需調(diào)用追趕法求解三對(duì)角線性方程組的算法 ? ? case 'CrankNicolson' ? ? ?? ?Diag=zeros(1,M-1);%矩陣的對(duì)角線元素 ? ? ?? ?Low=zeros(1,M-2);%矩陣的下對(duì)角線元素 ? ? ?? ?Up=zeros(1,M-2);%矩陣的上對(duì)角線元素 ? ? ?? ?for i=1:M-2 ? ? ?

35、? ?? ? Diag(i)=4; ? ? ?? ?? ? Low(i)=-r*C; ? ? ?? ?? ? Up(i)=r*C; ? ? ?? ?end ? ? ?? ?Diag(M-1)=4; ? ? ?? ?B=zeros(M-1,M-1); ? ? ?? ?for i=1:M-2 ? ? ?? ?? ? B(i,i)=4; ? ? ?? ?? ? B(i,i+1)=-r*C; ? ? ?? ?? ? B(i+1,i)=r*C; ? ? ?? ?end ? ? ?? ?B(M-1,M-1)=4; ? ? ?? ?%逐層求解，需要使用追趕法（調(diào)用函數(shù)EqtsForwardAndBack

36、ward） ? ? ?? ?for j=1:N ? ? ?? ?? ? b1=zeros(M-1,1); ? ? ?? ?? ? b1(1)=r*C*(U(1,j+1)+U(1,j))/2; ? ? ?? ?? ? b1(M-1)=-r*C*(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); ? ? ?? ?end ? ? ?? ? ? ? otherwise ? ? ?? ?disp('差分格式類型輸入有誤！') ? ? ?? ?return; end U=U'; %作出圖形 mesh(x,t,U); title([type '格式求解一階雙曲型方程的解的圖像']); xlabel('空間變量 x'); ylabel('時(shí)間變量 t'); zlabel('一階雙曲型方程的解 U'); return; 專心---專注---專業(yè)