使用 matlab 有限差分法在圆盘形域上的 Poisson PDE 求解器

Question

为了我的研究，我必须使用有限差分法为圆盘形域上的泊松方程编写 PDE 求解器。
我已经通过了实验室练习。我的代码中有一个问题我无法修复。具有边界值问题gun2 的函数fun1 在边界处以某种方式振荡。当我使用fun2 时，一切似乎都很好......

两个函数都在边界处使用gun2。有什么问题？

function z = fun1(x,y)
    r = sqrt(x.^2+y.^2);
    z = zeros(size(x));
    if( r < 0.25)
        z = -10^8*exp(1./(r.^2-1/16)); 
    end   
end 

function z = fun2(x,y)
    z = 100*sin(2*pi*x).*sin(2*pi*y);
end

function z = gun2(x,y)
   z = x.^2+y.^2;
end

function [u,A] = poisson2(funame,guname,M)

if nargin < 3
    M = 50;
end


%Mesh Grid Generation
h = 2/(M + 1);
x = -1:h:1;
y = -1:h:1;
[X,Y] = meshgrid(x,y);

CI = ((X.^2 +Y.^2) < 1);

%Boundary Elements
Sum= zeros(size(CI));

%Sum over the neighbours
for i = -1:1
    Sum = Sum + circshift(CI,[i,0]) + circshift(CI,[0,i]) ;
end

%if sum of neighbours larger 3 -> inner note!
CI = (Sum > 3);
%else boundary
CB = (Sum < 3 & Sum ~= 0);

Sum= zeros(size(CI));

%Sum over the boundary neighbour nodes....
for i = -1:1
   Sum = Sum + circshift(CB,[i,0]) + circshift(CB,[0,i]);
end

%If the sum is equal 2 -> Diagonal boundary
CB = CB + (Sum == 2 & CB == 0 & CI == 0);

%Converting X Y to polar coordinates
Phi = atan(Y./X);

%Converting Phi R back to cartesian coordinates, only at the boundarys
for j = 1:M+2
    for i = 1:M+2
        if (CB(i,j)~=0)
            if j > (M+2)/2
                sig = 1;
            else
                sig = -1;
            end
            X(i,j) = sig*1*cos(Phi(i,j));
            Y(i,j) = sig*1*sin(Phi(i,j));
        end
    end
end




 %Numberize the internal notes u1,u2,......,un

 CI = CI.*reshape(cumsum(CI(:)),size(CI));

 %Number of internal notes

 Ni = nnz(CI);


 f = zeros(Ni,1);

 k = 1;

 A = spalloc(Ni,Ni,5*Ni);



 %Create matix A!
 for j=2:M+1
    for i =2:M+1
        if(CI(i,j) ~= 0)
            hN = h;hS = h; hW = h; hE = h;

            f(k) = fun(X(i,j),Y(i,j));                                    
            if(CB(i+1,j) ~= 0)                                          
                hN = abs(1-sqrt(X(i,j)^2+Y(i,j)^2));
                f(k) = f(k) + gun(X(i,j),Y(i+1,j))*2/(hN^2+hN*h);
                A(k,CI(i-1,j)) = -2/(h^2+h*hN);
            else

                if(CB(i-1,j) ~= 0)  %in negative y is a boundry
                    hS = abs(1-sqrt(X(i,j)^2+Y(i,j)^2));
                    f(k) = f(k) + gun(X(i,j),Y(i-1,j))*2/(hS^2+h*hS);
                    A(k,CI(i+1,j)) = -2/(h^2+h*hS);
                else
                    A(k,CI(i-1,j)) = -1/h^2;
                    A(k,CI(i+1,j)) = -1/h^2;
                end
            end

            if(CB(i,j+1) ~= 0)                                          
                 hE = abs(1-sqrt(X(i,j)^2+Y(i,j)^2));
                 f(k) = f(k) + gun(X(i,j+1),Y(i,j))*2/(hE^2+hE*h);
                 A(k,CI(i,j-1)) = -2/(h^2+h*hE);
            else

                if(CB(i,j-1) ~= 0)
                     hW = abs(1-sqrt(X(i,j)^2+Y(i,j)^2));

                     f(k) = f(k) + gun(X(i,j-1),Y(i,j))*2/(hW^2+h*hW);
                     A(k,CI(i,j+1)) = -2/(h^2+h*hW);
                else
                     A(k,CI(i,j-1)) = -1/h^2;
                     A(k,CI(i,j+1)) = -1/h^2;
                end
            end

            A(k,k) = (2/(hE*hW)+2/(hN*hS));

            k = k + 1;
        end
    end
end

%Solve linear system
u = A\f;
U = zeros(M+2,M+2);
p = 1;

%re-arange u
for j = 1:M+2
    for i = 1:M+2
        if ( CI(i,j) ~= 0)
            U(i,j) = u(p);
            p = p+1;
        else
            if ( CB(i,j) ~= 0)
                U(i,j) = gun(X(i,j),Y(i,j));
            else
                U(i,j) = NaN;
            end
        end
    end
end


surf(X,Y,U);

end

Answer 1

我暂时保留这个答案，但当问题包含更多信息时可能会延长。

我的第一个猜测是，您看到的只是数字错误。查看两张图的比例，与第二张图中的信号相比，第一张图中的峰值相对较小。也许第二个中存在类似的问题，因为信号更大，所以不可见。您可以尝试增加节点数并观察结果。

您应该始终期望在此类模拟中看到数值错误。只是尝试使它们的大小尽可能小（或尽可能小）。