I wrote a Matlab code for solving 2D Poisson equation $u_{xx} + u_{yy} + f(x,y) = 0$ on $[a,b]\times [c,d]$ with neumann boundary condition on $x = b$ and the other boundary conditions are dirichlet,I used finite difference methods.My code is as follows:
function [u,err] = poissonFDM(a,b,c,d,j1,j2,fun,g,typebound,bound,uex,nit)
dltx = (b - a)/(j1+1); dlty = (d - c)/(j2+1);
switch typebound
case 'Neumann'
dim = (j1+1)*j2;
D = speye(j1+1);
T = ones(j1+1,1);
T = [T,-4*T,T]; T(j1,1) = 2*T(j1,1);
T = spdiags(T,[-1;0;1],j1+1,j1+1);
A1 = D;A = T;
for i = 1 :((j2)- 2)
A = blkdiag(A,T);
A1 = blkdiag(A1,D);
end
A = blkdiag(A,T);
A(1:(j1+1)*(j2-1),j1+2:(j1+1)*j2) = A(1:(j1+1)*(j2-1),j1+2:(j1+1)*j2) + A1;
A(j1+2:(j1+1)*j2,1:(j1+1)*(j2-1)) = A(j1+2:(j1+1)*j2,1:(j1+1)*(j2-1)) + A1;
h = max(dltx,dlty);
xh = (a:dltx:b); yh = (c:dlty:d)';
boundL = g(a,yh);
boundB = g(xh,c);boundU = g(xh,d);
gR = bound(b,yh);
k = 1;
v = zeros(dim,1);f = zeros(dim,1);Uex = zeros((j1+1)*j2,1);
for s = 1:j2
for r = 1:(j1+1)
f(r,s) = fun(r*dltx + a,s*dlty + c);
if(nargin == 12)
Uex(k) = uex(r*dltx+a,s*dlty+c);
end
v(k) = -(h^2)*f(r,s);
k = k + 1;
end
end
%Bottom Border
v(1:j1) = v(1:j1) - boundB(2:j1+1)';
v(j1+1) = v(j1+1) - gR(1);
%Left Border
v(1:j1:j1*j2 - j1 + 1) = v(1:j1:j1*j2 - j1 + 1) - boundL(2:j2+1);
%Right Border
v(j1+1:j1+1:(j1+1)*j2) = v(j1+1:j1+1:(j1+1)*j2) - 2*h*gR(2:j2+1);
%Up Border
v((j1+1)*j2-j1:(j1+1)*j2 -1) = v((j1+1)*j2-j1:(j1+1)*j2-1) - boundU(2:j1+1)';
v((j1+1)*j2)= v((j1+1)*j2 ) - gR(j2+2);
%solving linear system with conjugate gradiant
%U = conjgrad(A,v);
U = A\v;
%err = U - Uex;
end
end
I doubt about accuracy of my code ,I suppose that $u = \sin(x+y)$.when I test my code with this :
u = poissonFDM(0,pi,0,pi,11,11,@(x,y)2*sin(x+y),@(x,y)sin(x+y),'Neumann',@(x,y)cos(x+y),@(x,y)sin(x+y));
surf(u)
and this graph is far from the graph of $ u = sin(x+y)$. What is wrong with my code ??and about the method I used ,I should say that I approximate all the derivatives with the central finite difference formula!Can anyone help me?