偏微分方程数值解

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