《偏微分方程數(shù)值解法的MATLAB源碼》由會員分享，可在線閱讀，更多相關(guān)《偏微分方程數(shù)值解法的MATLAB源碼（16頁珍藏版）》請在裝配圖網(wǎng)上搜索。
1、精選優(yōu)質(zhì)文檔-----傾情為你奉上
[ 原創(chuàng)]偏微分方程數(shù)值解法的MATLAB源碼【更新完畢】
說明：由于偏微分的程序都比較長，比其他的算法稍復(fù)雜一些，所以另開一貼，專門上傳偏微分的程序
謝謝大家的支持！
其他的數(shù)值算法見：
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 -時間變量
%輸入?yún)?shù)：uX -空間變量x的取值上限
% ? ?? ?? uT -時間變量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
%計
4、算步長
dx=uX/M;%x的步長
dt=uT/N;%t的步長
x=(0:M)*dx;
t=(0:N)*dt;
r=C*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));
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('時間變量 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 -時間變量
%輸入?yún)?shù)：uX -空間變量x的取值上限
% ? ??
7、 ?? 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ù)
% ? ?? ?? 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
%計算步長
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
9、;
end
Diag(M-1)=1+2*r;
%計算初值和邊值
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('時間變量 t')
zlabel('一維熱傳導(dǎo)方程的解 U')
return;
此算法需要使用追趕法求解三對角線性方程組，這個算法在上一篇帖子中已經(jīng)給出，為了方便，再給出來
追趕法解三對角線性方程組
function x=EqtsForwardAndBackward(L,D,U,b)
%追趕法求解三對角線性方程組Ax=
11、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);
if n-n1 ~= 1
12、 || 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);
end
x=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(n);
for i=n-1:-1:1
? ? x(i)=(x(i)-U(i)*x(i+1))/D(i);
end
return;
13、
古典隱式格式
在以后的程序中，我們都取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 <= 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 -時間變量
%輸入?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ù)
% ? ?? ?? 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);
%計算步長
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
?
16、 ? Diag(i)=1+r;
? ? Low(i)=-r/2;
? ? Up(i)=-r/2;
end
Diag(M-1)=1+r;
%計算初值和邊值
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('時間變量 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ū)域Laplace方程的Diriclet邊值問題的差分求解
%此程序需要調(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時的值，第二行表示第y=h時的值……
% ? ?? ?? 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
%步長
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('差分格式類型輸入錯誤')
? ? ?? ?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問題的解的圖像')
xlabel('x')
ylabel('y')
zlabel('Laplace
25、方程Diriclet問題的解 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個時間層……
% ? ?? ? x -橫坐標(biāo)
% ? ?? ? t -縱坐標(biāo)，時間
%輸入?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)用追趕法
% ? ?? ?? ?? ? 求解三對角線性方程組
%
h=uX/M;%變量x的步長
k=uT/N;%變量t的步長
r=
28、k/h;%步長比
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('請輸入第二層初值條件函數(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)用追趕法求解三對角線性方程組的算法
? ? case 'CrankNicolson'
? ? ?? ?Diag=zeros(1,M-1);%矩陣的對角線元素
? ? ?? ?Low=zeros(1,M-2);%矩陣的下對角線元素
? ? ?? ?Up=zeros(1,M-2);%矩陣的上對角線元素
? ? ?? ?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('時間變量 t');
zlabel('一階雙曲型方程的解 U');
return;
專心---專注---專業(yè)
偏微分方程數(shù)值解法的MATLAB源碼
