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

[omer21]

I am trying to solve the following pde numerically using backward f.d. for time and central difference approximation for x, in MATLAB but i can't get correct results.

$\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$
for boundary conditions i used the following approximation
$u_{x}(0,t)=\frac{u_{i+1}^{j}-u_{i-1}^{j}}{2h}$

what is wrong with the code i wrote

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,:);

 

[jfgobin]

Hey Omer21,

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

$u\left ( 0,0 \right ) =2\\ u\left ( 1,0\right ) = 1$
​

J.

 

[omer21]

Initial vector is right. You miss subscript at b. c. i guess. B.C's are specified at the derivative of u.

 

[jfgobin]

Whatever. If you want a discontinuous function as boundary condition, that's your right.

Without having looked further, I see you are missing a

T(m+1,j)=2;

in your main loop. Also, the following lines are suspicious:

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.

 

[omer21]

You are right about

T=T(:,1);

. It ought to be

T=T(:,end);

.
I want to show what i did so i multiplied by 0 because of $u'_x(0,t)=0$.
I put

T(m+1,j)=2;

just after second for loop but still i can't get correct results.

 

[jfgobin]

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.

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

 

[omer21]

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.

 

[jfgobin]

Okay, I think I know where the problem lies. Part here, part there.

I will write some code later.

 

[jfgobin]

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.

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?

 

[omer21]

Your plots are correct. Thank you for your effort.