1. Limited time only! Sign up for a free 30min personal 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!

Homework Help: 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:


    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
    d1=-ones(N*N,1);            %sub and superdiagonals

    for i=N:N:N*N-1

    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
    for i=1:N                   %inner points as a matrix
    U=zeros(N+2,N+2);           %add BCs to inner solution
    for i = 1:N
    mesh(x,y,U)                 %3D plot of the solution
    title('Solution of \Laplace u = 1 on the unit square')
    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

    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
    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.

    Last edited: Nov 11, 2013
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted