0
$\begingroup$

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 : enter image description here

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?

$\endgroup$
  • 2
    $\begingroup$ It might help if you explain (in equations) how you implemented your scheme and your boundary conditions. Then explain how you translated the math into your code. $\endgroup$ – Paul Jul 1 '16 at 16:42
1
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
0
$\begingroup$

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 !

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Never mind 3, you commented 'Uex' out on purpose. My bad. :) $\endgroup$ – Inquisitor101 Jul 6 '16 at 11:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.