偏微分方程數值解

1.求解拉普拉斯方程的狄利克雷法

求解在區域R = {(x,y): 0≤x≤a, 0≤y≤b}內的 u_xx(x,y) + u_yy(x,y) = 0 的近似解，並且知足條件 u(x,0) = f₁(x),  u(x,b) = f₂(x), 其中0≤x≤a 且 u(0,y) = f₃(y), u(a,y) = f₄(y),其中 0≤y≤b。設Δx = Δy = h，並且存在整數n和m，使得 a = nh ,b = mh。

代碼以下：

function [U,cnt]=dirich(f1,f2,f3,f4,a,b,h,tol,max1)
%Input  - f1,f2,f3,f4 are the function entered  as a string 
%       - a and b are the left and right end points
%       - h steps size
%       - tol is the tolerance
%       - max1 is maximum of iterations number
%Output - U solution matrix;analogous to Table 10.6
%         cnt is number of iterations
%
%其算法的思路是，先將內部網格點均取四個角點的平均值做爲內部全部網格點的初值，
%而後用SQR超鬆弛算法
%內部網格點的迭代值爲初值點+餘項，餘項爲拉普拉斯差分計算方程
%對每一個網格點不斷迭代，使得最終餘項 relx趨近於0爲止（知足拉普拉斯差分方程等於零特色）。
%Initialize parameters and U
n=fix(a/h)+1;
m=fix(b/h)+1;
ave=(a*(feval(f1,0)+feval(f2,0))+b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b);
U=ave*ones(n,m);

%Boundary conditions
U(1,1:m)=feval(f3,0:h:(m-1)*h)';
U(n,1:m)=feval(f4,0:h:(m-1)*h)';
U(1:n,1)=feval(f1,0:h:(n-1)*h);
U(1:n,m)=feval(f2,0:h:(n-1)*h);

%SOR parameter
w=4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));

%Refine approximations and sweep operator throughout the grid
err=1;
cnt=0;
while((err>tol)&&(cnt<=max1))
    err=0;
    for j=2:m-1
        for i=2:n-1
            relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j))/4;
            U(i,j)=U(i,j)+relx;
            if (err<=abs(relx))
                err=abs(relx);
            end
        end
    end
    cnt=cnt+1;
end
U=flipud(U');%flipud實現矩陣的上下翻轉
end