Poisson偏微分方程求解器使用matlab有限差分法求解形狀區域

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

我保持這樣的回答簡短現在，但是當問題包含更多信息可擴展。

我的第一個猜測是，你所看到的只是數字錯誤。查看兩幅圖的比例尺，第一幅圖中的峯與第二幅圖中的信號相比相對較小。也許第二個類似的問題只是不可見，因爲信號更大。您可以嘗試增加節點的數量並觀察結果。

您應該始終期望在這種模擬中看到數值錯誤。這只是一個試圖儘可能小（或根據需要儘可能小）的問題。