此處參考教材為李榮華的《微分方程數值解法》
使用工具：Matlab

1. 算法：矩形網格上5點差分格式

2. 算法

I.需要求解的函數

function [v,vx,vy,f,aa,bb,cc,dd]=u2D(x,y,ft)% ft為方程編號，u1D為精確解函數u（t）,注意與f2D 對應右端項函數f(t,u(t))switch ftcase 1 %P-104例題v=cos(3*x).*sin(pi*y)/(9+pi^2);vx=-3*sin(3*x).*sin(pi*y)/(9+pi^2);vy=pi*cos(3*x).*cos(pi*y)/(9+pi^2);f=cos(3*x).*sin(pi*y);aa=0;bb=pi;cc=0;dd=1;case 2 %P-104習題v=exp(pi*(x+y)).*sin(pi*x).*sin(pi*y); vx=pi*exp(pi*(x+y)).*(sin(pi*x).*sin(pi.*y)+cos(pi*x).*sin(pi*y)); vy=pi*exp(pi*(x+y)).*(sin(pi*x).*sin(pi*y)+cos(pi*y).*sin(pi*x)); f=-2*(pi^2)*exp(pi*(x+y)).*(sin(pi*x).*cos(pi*y)+cos(pi*x).*sin(pi*y)); aa=0;bb=1;cc=0;dd=1; end

II.求解例題的主函數：五點差分格式函數

function Elliptic_PDESolver(ft,M,N) %初始化定義求解區域和網格 [~,~,~,~,x1,xr,y1,yr]=u2D(0,0,ft); %求解區域也從u2D函數調用，初始x和y設為0 h1=(xr-x1)/M; h2=(yr-y1)/N; xd=(x1:h1:xr); %網格節點坐標 yd=(y1:h2:yr); %網格節點坐標%定義全部節點編號集AId，網格節點編號與矩陣下標編號不一致 AId=zeros(M+1,N+1); xk=zeros((M+1)*(N+1),1); %整體編號排序后的節點坐標 yk=zeros((M+1)*(N+1),1); %整體編號排序后的節點坐標 for i=0:Mfor j=0:NAId(i+1,j+1)=i*(N+1)+j+1; %注意矩陣編號從1開始，故i+1，j+1xk(i*(N+1)+j+1)=xd(i+1); %xd的編號從1開始，故i+1yk(i*(N+1)+j+1)=yd(j+1); %yd的編號從1開始，故j+1end end AId=AId';AId=AId(:);%定義內部節點編號集：InId InId=zeros(M-1,N-1); for i=1:M-1for j=1:N-1InId(i,j)=i*(N+1)+j+1;end end InId=InId';InId=InId(:);%定義左邊界節點編號集 LbId for j=0:NLbId(j+1)=0*(N+1)+j+1; end%定義右邊界節點編號集 RbId for j=0:NRbId(j+1)=M*(N+1)+j+1; end%定義下邊界節點編號集 　BbId for i=0:MBbId(i+1)=i*(N+1)+0+1; end%定義上邊界節點編號集 　TbId for i=0:MTbId(i+1)=i*(N+1)+N+1; end%顯示整體節點編號圖 figure; for i=0:Mplot(xd(i+1)+0*yd,yd);hold on; end for j=0:Nplot(xd,yd(j+1)+0*xd);hold on; end for i=1:length(AId)text(xk(i),yk(i),num2str(i)); end%初始化矩陣右端項 A=zeros((M+1)*(N+1),(M+1)*(N+1)); b=zeros((M+1)*(N+1),1); hh1=h1^2; hh2=h2^2;%內部節點形成矩陣行元素及右端項 for i=1:length(InId)k=InId(i);A(k,k)=2*hh1+2*hh2;A(k,k-1)=-hh1;A(k,k+1)=-hh1; A(k,k-N-1)=-hh2;A(k,k+N+1)=-hh2;[~,~,~,ff]=u2D(xk(k),yk(k),ft); %提取右端項函數值b(k)=hh1*hh2*ff; end %添加邊界條件元素及右端項 %定義左邊界節點矩陣行元素及右端項 for j=1:length(LbId) k=LbId(j); if ft==1A(k,k)=-1;A(k,k+N+1)=1; [~,ux,~,~]=u2D(xk(k),yk(k),ft); %%提取右端項函數值，第二邊值條件b(k)=h1*ux; elseif ft==2A(k,k)=1;A(k,k+N+1)=0;[uv,~,~,~]=u2D(xk(k),yk(k),ft); %提取右端項函數值;第一邊值條件b(k)=uv;%b(k)=0;end end% 定義右邊界節點矩陣行元素及右端項 for j=1:length(RbId) k=RbId(j);if ft==1A(k,k)=-1;A(k,k-N-1)=1; [~,ux,~,~]=u2D(xk(k),yk(k),ft); %提取右端項函數數值，第二邊值條件b(k)=h1*ux;elseif ft==2A(k,k)=1;A(k,k-N-1)=0;[uv,~,~,~]=u2D(xk(k),yk(k),ft); %提取右端項函數數值，第一邊值條件b(k)=uv;%b(k)=0;end end% 定義下邊界節點矩陣行元素及右端項 for i=1:length(BbId)k=BbId(i);A(k,:)=0; %由于與其他邊界條件有重復，清空該行元素A(k,k)=1;[uv,~,~,~]=u2D(xk(k),yk(k),ft); %提取右端項函數值，第一邊值條件b(k)=uv;%b(k)=0; end% 定義上邊界節點矩陣行元素及右端項 for i=1:length(TbId)k=TbId(i);A(k,:)=0; %由于與其他邊界條件有重復，清空該行元素A(k,k)=1;[uv,~,~,~]=u2D(xk(k),yk(k),1); %提取右端項函數值，第一邊值條件b(k)=uv;%b(k)=0; enduk=A\b; %求解所得整體節點排序的數值解 ExUk=u2D(xk,yk,ft); %所得整體節點排序的精確解 MaxErr=max(abs(uk-ExUk)); fprintf('MaxErr=%3.6f\n',MaxErr);%將數值解轉化為網格排序 NuU=zeros(M+1,N+1); for i=0:Mfor j=0:Nk=i*(N+1)+j+1;NuU(i+1,j+1)=uk(k);end end NuU=NuU'; %轉置后矩陣下標與U相同%精確解畫圖 [X,Y]=meshgrid(xd,yd); U=u2D(X,Y,ft); figure; mesh(X,Y,U); xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis'); hold on; %數值解畫圖 plot3(xk,yk,uk,'+'); hold on;

III.補充算例：八邊形區域剖分

function Elliptic_PDESolver2(ft,M,N) fprintf('八邊形區域網格剖分數M和N需相等并且為3的倍數 \n'); if(M~=N)||(mod(M,3)~=0)fprintf('輸入M和N不正確 \n');return end%clear all % 初始化定義求解區域和網格 %ft=1; %M=18;%N=18;[~,~,~,~,x1,xr,y1,yr]=u2D(0,0,ft); %求解區域也從u2D函數調用，初始x和y設為0 h1=(xr-x1)/M; h2=(yr-y1)/N; xd=(x1:h1:xr); %網格節點坐標 yd=(y1:h2:yr); %網格節點坐標% 定義全部節點編號集 AId，注意網格節點編號與矩陣下標編號 AId=zeros(M+1,N+1); xk=zeros((M+1)*(N+1),1); %整體排序后的節點坐標 yk=zeros((M+1)*(N+1),1); %整體排序后的節點坐標 for i=0:Mfor j=0:NAId(i+1,j+1)=i*(N+1)+j+1; %注意矩陣編號從1開始，故i+1，j+1xk(i*(N+1)+j+1)=xd(i+1); %xd的編號從1開始，故i+1yk(i*(N+1)+j+1)=yd(j+1); %yd的編號從1開始，故i+1end end AId=AId';AId=AId(:);%顯示整體節點編號圖 figure; for i=0:Mplot(xd(i+1)+0*yd,yd,'-.');hold on; end for j=0:Nplot(xd,yd(j+1)+0*xd,'-.');hold on; end for i=1:length(AId)text(xk(i),yk(i),num2str(i)); end%八邊形區域 %定義內部節點編號集 InId InId=[]; B1Id=[]; %四條斜邊 B2Id=[]; %左邊 B3Id=[]; %右邊 B4Id=[]; %下邊 B5Id=[]; %上邊 for i=0:Mfor j=0:Nif(i>0)&&(i<M)&&(j>0)&&(j<N)&&(j>N/3-i)&&(j<2*N/3+i)&&(j>i-2*N/3)&&(j<5*N/3-i)InId=[InId,i*(N+1)+j+1];endif(j==N/3-i)||(j==2*N/3+i)||(j==i-2*N/3)||(j==5*N/3-i)B1Id=[B1Id,i*(N+1)+j+1];endif(i==0)&&(j>N/3)&&(j<2*N/3)B2Id=[B2Id,i*(N+1)+j+1];endif(i==M)&&(j>N/3)&&(j<2*N/3)B3Id=[B3Id,i*(N+1)+j+1];endif(j==0)&&(i>M/3)&&(i<2*M/3)B4Id=[B4Id,i*(N+1)+j+1];endif(j==N)&&(i>M/3)&&(i<2*M/3)B5Id=[B5Id,i*(N+1)+j+1];endend end plot(xk(InId),yk(InId),'ro'); hold on; plot(xk(B1Id),yk(B1Id),'b*'); hold on; plot(xk(B2Id),yk(B2Id),'g*',xk(B3Id),yk(B3Id),'g*'); hold on; plot(xk(B4Id),yk(B4Id),'c*',xk(B5Id),yk(B5Id),'c*'); hold on;% 初始化矩陣和右端項 A=zeros((M+1)*(N+1),(M+1)*(N+1)); b=zeros((M+1)*(N+1),1); hh1=h1^2; hh2=h2^2;%內部節點形成矩陣行元素及右端項 for i=1:length(InId)k=InId(i); %提取內部節點編號A(k,k)=2*hh1+2*hh2;A(k,k-1)=-hh1;A(k,k+1)=-hh1; A(k,k-N-1)=-hh2;A(k,k+N+1)=-hh2;[~,~,~,ff]=u2D(xk(k),yk(k),1); %提取右端項函數數值b(k)=hh1*hh2*ff;%b(k)=hh1*hh2*(cos(3*xk(k))*sin(pi*yk(k))); end%添加邊界條件矩陣元素及右端項 % 定義B1Id邊界節點矩陣行元素及右端項 for j=1:length(B1Id)k=B1Id(j);A(k,k)=1;[uv,~,~,~]=u2D(xk(k),yk(k),1); %提取右端項函數數值,第一邊值條件b(k)=uv; end% 定義B2Id邊界節點矩陣行元素及右端項 for j=1:length(B2Id)k=B2Id(j);A(k,k)=-1;A(k,k+N+1)=1;[~,ux,~,~]=u2D(xk(k),yk(k),1); %提取右端項函數數值，第二邊值條件b(k)=h1*ux; end% 定義B3Id邊界節點矩陣行元素及右端項 for j=1:length(B3Id)k=B3Id(j);A(k,k)=1;A(k,k-N-1)=-1;[~,ux,~,~]=u2D(xk(k),yk(k),1); %提取右端項函數數值，第二邊值條件b(k)=h1*ux; end% 定義B4Id邊界節點矩陣行元素及右端項 for j=1:length(B4Id)k=B4Id(j);A(k,k)=-1;A(k,k+1)=1;[~,~,uy,~]=u2D(xk(k),yk(k),1); %提取右端項函數數值，第一邊值條件b(k)=h2*uy; end% 定義B5Id邊界節點矩陣行元素及右端項 for j=1:length(B5Id)k=B5Id(j);A(k,k)=1;A(k,k-1)=-1;[~,~,uy,~]=u2D(xk(k),yk(k),1); %提取右端項函數數值，第一邊值條件b(k)=h2*uy; end%矩陣A是按照所有(M+1)*(N+1)個節點生成的矩陣，里面有很多的非區域節點編號對應的0子陣 DmId=sort([InId,B1Id,B2Id,B3Id,B4Id,B5Id]); OutId=setdiff(AId,DmId); A(OutId,:)=[];A(:,OutId)=[]; %剔除非區域節點OutId對應的行和列 b(OutId,:)=[]; %剔除非區域節點OutId對應的右端項 uk=A\b; %求解所得整體節點排序的數值解 ExUk=u2D(xk(DmId),yk(DmId),ft); %整體節點排序的精確解 MaxErr=max(abs(uk-ExUk)); fprintf('MaxErr=%3.6f\n',MaxErr);%精確解畫圖 [X,Y]=meshgrid(xd,yd); U=u2D(X,Y,ft); figure; mesh(X,Y,U); xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis'); hold on; %數值解畫圖 plot3(xk(DmId),yk(DmId),uk,'+'); hold on;

3. 給出簡單結果

[^此處圖略，僅給出簡單數字結果；感興趣的讀者可以復制代碼運行一下]

P104頁：習題

（1） 步長k=h=1/64：

在命令行窗口輸入
>>Elliptic_PDESolver(2,64,64)

所得結果：

MaxErr=0.028613

（2）步長k=h=1/128：

在命令行窗口輸入
>> Elliptic_PDESolver(2,128,128)
所得結果：

MaxErr=0.007156

P104頁：例題

（1） N=M=4：在命令行窗口輸入

>> Elliptic_PDESolver(1,4,4)
所得結果：

MaxErr=0.117314

之后重復上面命令行的步驟

（2）M=N=8

MaxErr=0.047304

（3）M=N=16

MaxErr=0.019118

（4）M=N=32

MaxErr=0.008457

補充練習

>> Elliptic_PDESolver2(1,12,12)

八邊形區域網格剖分數M和N需相等并且為3的倍數
MaxErr=0.019358

>> Elliptic_PDESolver2(1,15,15)

八邊形區域網格剖分數M和N需相等并且為3的倍數
MaxErr=0.013485

4. 對算例結果進行簡單的分析

在這兩道數值例子中，顯然當網格加密的時候誤差越來越小，精確程度越來越高。但是在matlab程序運行過程中明顯能感覺到，網格越密運行時間越長。
由補充練習知 八邊形區域網格剖分數M和N需相等并且為3的倍數，且隨著網格的加密誤差逐漸變小。

總結

以上是生活随笔為你收集整理的微分方程数值解法（2）——椭圆型方程的有限差分法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。