What is wrong with my code for solving Poisson equation with one side Neumann boundary condition?

Question

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)

so the graph will be this :

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?

Answer 1

It's unclear what the problem is to begin with. Why do you doubt the accuracy of your solution? If you substitute the numerical result into the governing equation do you return a zero residual? A central difference implementation should be fine for this problem, except at the boundaries of course.

Try simplifying your problem to find the error: remove the forcing term $f(x,y)$ so you're solving a Laplace equation, simplify the boundary conditions (Dirichlet on all sides perhaps), reduce the order so you're solving the 1D problem, and so fourth.

Answer 2

The question you wrote down is very puzzling. To start with, taking the same code and pasting it while using the same command you used to generate the output yields an error.

1 - You specified 'u' (lower-case) as a function output, yet it is used nowhere in your code (you only use 'U' -- upper case).

2 - When I substitute 'u' for 'U' -- as one expects -- I get a vector. Meaning I cannot reproduce the output via 'surf' as you did. Please synchronize your question for us to understand how to assist.

3 - Out of curiosity, 'Uex' is used where exactly after calculating it in the inner-most loop?

4 - Some other variables seem not be used at all, however I will assume this was done on purpose to simplify the question...

5 - Would it not be easier if you used a more simplified function -- taking into consideration how you'd expect the boundaries to react -- for debugging purposes ?

6 - Why are there no comments in your code. This is a bad programming habit, even for a relatively simple function and not to mention an inquiry post.

To sum it up, consider monitoring the residual (every $k^{th}$ iteration?), simplify your input function and reduce your boundaries to Dirichlet. An even more useful tool is to insert a breakpoint(s) inside your function, to further break down your code and debug effectively. If this still would not work, consider fixing your question by making it reproducible on our side and please: comments

Good luck !