1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Using finite difference method for solving an elliptic PDE with MATLAB

  1. Nov 11, 2013 #1
    1. The problem statement, all variables and given/known data

    Given that we the following elliptic problem on a rectangular region:

    [itex]\nabla^2 T=0, \ (x,y)\in \Omega[/itex]

    [itex]T(0,y)=300, \ T(4,y)=600, \ 0 \leq y \leq 2[/itex]

    [itex]\frac{\partial T}{\partial y}(x,0)=0, \frac{\partial T}{\partial y}(x,2) = 0, \ 0\leq x \leq 4[/itex]

    We want to solve this problem using finite difference method with stepsize h = 0.2, M divisions on x-axis and N on y-axis.
    My question is how can i find the Neuman BC in "discretized form?
    Since we need to exclude the boundaries, how many internal unknown equations will there be inside the region?
    Also, exactly how should the implementation be in Matlab code?


    2. Relevant equations

    I found the discretizion of the Laplace equation to be (using standard difference formulas for the equation), after simplification:

    [itex]-4T_{i,j}+T_{i-1,j}+T_{i+1,j}+T_{i,j-1}+T_{i,j+1}=0[/itex]

    For the matlab code, I have an example code that solves [itex]\nabla^2 u = 1[/itex] with BC [itex]u=0[/itex] at all boundaries.

    Code (Text):
    clear all;

    N = 100;                    %Number of inner x and y-points
                               
    hx=1/(N+1); hy=1/(N+1);     %stepsize in x- and y-direction

    x=0:hx:1; y=0:hy:1;         %gridpoints
    d=4*ones(N*N,1);            %main diagonal of A
    d_1=-ones(N*N,1);
    d1=-ones(N*N,1);            %sub and superdiagonals

    for i=N:N:N*N-1
        d_1(i)=0
        d1(i+1)=0
    end

    d_2=-ones(N*N,1);           %subdiagonal from unit matrix
    d2=-ones(N*N,1);            %superdiagonal
    a = spdiags([d_2 d_1 d d1 d2], [-N, -1:1,N],N*N,N*N);  
                                %sparse matrix A generated
    f=ones(N*N,1);              %right hand side
    uinner=a\f;                 %solution at inner points
    Uinner=zeros(N,N);          
               
    for i=1:N                   %inner points as a matrix
        Uinner(i,1:N)=uinner((i-1)*N+[1:N]);
    end
    U=zeros(N+2,N+2);           %add BCs to inner solution
    for i = 1:N
        U(i+1,2:N+1)=Uinner(i,1:N);
    end
    mesh(x,y,U)                 %3D plot of the solution
    title('Solution of \Laplace u = 1 on the unit square')
        xlabel('x')
        ylabel('y')
        zlabel('u(x,y)')
    3. The attempt at a solution

    Since in the problem we cannot include the boundary points, I calculate, the inner points as [itex](N-1)(M-1)[/itex] which is all the unknown equations to be solved right?

    For the Dirichlet BC: i set that:

    [itex]T_{0,j}=300, \ j=1,2,3...N[/itex] and
    [itex]T_{4,j}=600, \ j=1,2,3...N[/itex]

    Now for the Neuman ones i just by analogy from lecture notes write down the conditions, but need to get confirmation if it is correct that:

    [itex]\frac{T_{i,0}-T_{i,2}}{2h}=0[/itex] and
    [itex]\frac{T_{i,0}-T_{i,2}}{2h}=0[/itex]

    Can anyone confirm if neuman conditions are done right?

    Also, in the matlab code I know that:
    Code (Text):
    for i=N:N:N*N-1
        d_1(i)=0
        d1(i+1)=0
    end
    corresponds to the boundary conditions but I do not understand what this for-loop actually does because at the end the -ones is still a vector with N*N rows with -1

    It was a long question..
    Hope someone can help me.

    Thanks
     
    Last edited: Nov 11, 2013
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?



Similar Discussions: Using finite difference method for solving an elliptic PDE with MATLAB
Loading...