Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Backward euler method for heat equation with neumann b.c.

  1. Mar 17, 2013 #1
    I am trying to solve the following pde numerically using backward f.d. for time and central di fference approximation for x, in matlab but i can't get correct results.

    [itex]
    \frac{\partial u}{\partial t}=\alpha\frac{\partial^{2}u}{\partial x^{2}},\qquad u(x,0)=f(x),\qquad u_{x}(0,t)=0,\qquad u_{x}(1,t)=2
    [/itex]
    for boundary conditions i used the following approximation
    [itex]u_{x}(0,t)=\frac{u_{i+1}^{j}-u_{i-1}^{j}}{2h}[/itex]

    what is wrong with the code i wrote
    Code (Text):

    function [T,exact]=implicitheat2(t_i,t_f,a,b,dx,dt,alpha)

    %U_t=U_xx , 0<x<1, U(x,0)=x.^2+1+cos(pi.*x), U_x(0,t)=0, U_x(1,t)=2

    % dt: step size in t
    % dx: step size in x
    % a: left point of domain
    % b: right point of domain
    % alpha: equal to 1
    % call func. as implicitheat2(0,0.1,0,1,0.1,0.00001,1)

     

    n=(t_f-t_i)/dt;
    m=(b-a)/dx;
    lambda=alpha*dt/dx^2;

    if isinteger(m)==0
        m=round(m);
    end
    if isinteger(n)==0
        n=round(n);
    end
    T=zeros(m+1,n+1);

    x=a:dx:b;
    t=t_i:dt:t_f;
       
    u0 = x.^2+1+cos(pi.*x);
    T(:,1)=u0; %initial value

    A = sparse(m-1,m-1);
        for i=1:m-1
            A(i,i-1) = -lambda;
            A(i,i ) = (1+2*lambda);
            A(i,i+1) = -lambda;
        end
    A(1,2)=-2*lambda;
    A(end,end-1)=-2*lambda;
    b=zeros(m-1,1);
    for j=2:n+1
        b=T(2:m,j-1);
        b(1,1)=T(2,j-1)-2*0*lambda*dx;  %T(0,j-1)-2*lambda*U'*dx
        b(m-1,1)=T(m,j-1)+2*lambda*2*dx;  %T(m+1,j-1)+2*lambda*U'*dx
        T(2:m,j)=A\b;
    end
    T=T(:,1);
    %exact soln
    [xx,tt]=meshgrid(x,t);

    exact=2.*tt+xx.^2+1+exp(-pi^2.*tt).*cos(pi.*xx);

    exact=exact(end,:);
     
     
  2. jcsd
  3. Mar 18, 2013 #2
    Hey Omer21,

    I will have a look at this a bit later. Also, could you check your initial vector [itex]u_0[/itex]? I plotted it and I have

    [itex]
    u\left ( 0,0 \right ) =2\\
    u\left ( 1,0\right ) = 1
    [/itex]
    J.
     
  4. Mar 18, 2013 #3
    Initial vector is right. You miss subscript at b. c. i guess. B.C's are specified at the derivative of u.
     
  5. Mar 18, 2013 #4
    Whatever. If you want a discontinuous function as boundary condition, that's your right.

    Without having looked further, I see you are missing a

    Code (Text):

    T(m+1,j)=2;
     
    in your main loop. Also, the following lines are suspicious:

    Code (Text):


    b(1,1)=T(2,j-1)-2*0*lambda*dx;

    T=T(:,1);
     
    First: you multiply by 0. You sure? And for the second, you are doing a bunch of calculations and then taking only the first column.
     
  6. Mar 18, 2013 #5
    You are right about
    Code (Text):
    T=T(:,1);
    . It ought to be
    Code (Text):
    T=T(:,end);
    .
    I want to show what i did so i multiplied by 0 because of [itex]u'_x(0,t)=0[/itex].
    I put
    Code (Text):
    T(m+1,j)=2;
    just after second for loop but still i can't get correct results.
     
  7. Mar 18, 2013 #6
    I redid your code. I think there is more to it than "just" backward Euler not working.

    As I mentioned -

    - Your initial vector first and last elements are not 0 and 2;
    - You are assigning weird things in weird ways;
    - Your "exact" solutions doesn't have value 0 for x=0 and value 2 for x=1, which are your space boundaries.

    Here is my code in Octave. Feel free to inspire yourself. I will add the averaging with forward Euler to increase the precision.

    Code (Text):

    function Y=heattrans(t0,tf,n,m,alpha,withfe)

    # Calculate the heat distribution along the domain 0->1 at time tf, knowing the initial
    # conditions at time t0
    # n - number of points in the time domain (at least 3)
    # m - number of points in the space domain (at least 3)
    # alpha - heat coefficient
    # withfe - average backward Euler and forward Euler to reach second order

    # The equation is
    #
    #   du             d2u
    #  ---- = alpha . -----
    #   dt             dx2

    if (or(n<3,m<3))
     disp("Really mate?!")
     Y=[]
     return
    endif

    dx=1/(m-1);
    dt=(tf-t0)/(n-1);

    # Initial conditions are
    # U(x,0) = x.^2 + 1 + cos(pi*x)
    # Warning - U(0,0) = 2 =/= U(0,t=/=0) = 0
    # and     - U(1,0) = 3 =/= U(1,t=/=0) = 2

    x=linspace(0,1,m);
    t=linspace(t0,tf,n);
    beta=(alpha*dt)/(dx^2);

    # Let's initialize our values
    Ut=x.^2+1+cos(pi*x);

    # Main loop - We apply backward Euler and solve successive equations
    # Second order differences are computed using central difference, backward Euler is computed
    # using the classic "one side" difference

    for k = 1:n
      # Let's build the matrix of diff factors
      M=spalloc(m,m,3*m);
      M(1,1)=1;
      M(m,m)=1;
      for l = 2:(m-1)
        M(l,l-1)=beta;
        M(l,l)=-(2*beta+1);
        M(l,l+1)=beta;
      endfor
      # Now, let's build the vector of diff terms
      # That's actually -Ut, except for the first and last elements
      D=-Ut';
      D(1,1)=0;
      D(m,1)=2;
      # And the next step is the solution (transposed, as I use line vectors for U)
      Ut=(inv(M)*D)';
    endfor
    # The answer is in the last Ut
    Y=Ut;
    endfunction
     
     
  8. Mar 18, 2013 #7
    Actually my codes are working accurately when B.C.s are dirichlet, but when B.C.s are turned to Neumann i confused how to edit the codes anyway i will look your codes.
    Thank you jfgobin for your help.
     
  9. Mar 19, 2013 #8
    Okay, I think I know where the problem lies. Part here, part there.

    I will write some code later.
     
  10. Mar 19, 2013 #9
    Here is the code that solves it. Don't forget that both backward Euler and forward Euler are methods of the first order, and that imprecision can creep up.

    My code doesn't use central difference for the first order derivative: the only cases I need them is for the corners. A better approximation could be made by taking more points.

    Code (Text):

    function Y=heattrans(t0,tf,n,m,alpha,withfe)

    # Calculate the heat distribution along the domain 0->1 at time tf, knowing the initial
    # conditions at time t0
    # n - number of points in the time domain (at least 3)
    # m - number of points in the space domain (at least 3)
    # alpha - heat coefficient
    # withfe - average backward Euler and forward Euler to reach second order

    # The equation is
    #
    #   du             d2u
    #  ---- = alpha . -----
    #   dt             dx2

    if (or(n<3,m<3))
     disp("Really mate?!")
     Y=[]
     return
    endif

    dx=1/(m-1);
    dt=(tf-t0)/(n-1);

    # Initial conditions are
    # U(x,0) = x.^2 + 1 + cos(pi*x)
    # Warning - U(0,0) = 2 =/= U(0,t=/=0) = 0
    # and     - U(1,0) = 3 =/= U(1,t=/=0) = 2

    x=linspace(0,1,m);
    t=linspace(t0,tf,n);
    beta=(alpha*dt)/(dx^2);

    # Let's initialize our values
    Ut=x.^2+1+cos(pi*x);

    # Main loop - We apply backward Euler and solve successive equations
    # Second order differences are computed using central difference, backward Euler is computed
    # using the classic "one side" difference

    for k = 1:n
      # Let's build the matrix of diff factors
      M=spalloc(m,m,3*m);
      M(1,1)=-1;
      M(1,2)=1;
      M(m,m)=1;
      M(m,m-1)=-1;
      for l = 2:(m-1)
        M(l,l-1)=beta;
        M(l,l)=-(2*beta+1);
        M(l,l+1)=beta;
      endfor
      # Now, let's build the vector of diff terms
      # That's actually -Ut, except for the first and last elements
      D=-Ut';
      D(1,1)=0;
      D(m,1)=2*dx;
      # And the next step is the solution (transposed, as I use line vectors for U)
      Ut=(inv(M)*D)';
    endfor
    # The answer is in the last Ut
    Y=Ut;
    endfunction
     
    Attached are two plots, one at t=0.05 and one at t=1. Is that what you are looking for?
     

    Attached Files:

  11. Mar 19, 2013 #10
    Your plots are correct. Thank you for your effort.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Backward euler method for heat equation with neumann b.c.
Loading...