title: Poisson PDE solver on a disked shaped domain with finite difference method using matlab
body: <p>For my studies I had to write a PDE solver for the Poisson equation on a disc shaped domain using the finite difference method.<br>
  I already passed the Lab exercise. There is one issue in my code I couldn't fix. Function <code>fun1</code> with the boundary value problem <code>gun2</code> is somehow oscillating at the boundary. When I use <code>fun2</code> everything seems fine...   </p>
  
  <p>Both functions use at the boundary <code>gun2</code>. What is the problem?</p>
  
  <p><a href="https://i.stack.imgur.com/DCwNB.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DCwNB.jpg" alt="function 1 oscillating at the boundary"></a>
  <a href="https://i.stack.imgur.com/DBuOX.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DBuOX.jpg" alt="function 2 no oscillating"></a></p>
  
  <pre><code>function z = fun1(x,y)
      r = sqrt(x.^2+y.^2);
      z = zeros(size(x));
      if( r &lt; 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 &lt; 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) &lt; 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 -&gt; inner note!
  CI = (Sum &gt; 3);
  %else boundary
  CB = (Sum &lt; 3 &amp; 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 -&gt; Diagonal boundary
  CB = CB + (Sum == 2 &amp; CB == 0 &amp; 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 &gt; (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
  </code></pre>