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

1. Mar 17, 2013

### omer21

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.

$\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
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. Mar 18, 2013

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

3. Mar 18, 2013

### omer21

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

4. Mar 18, 2013

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

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.

5. Mar 18, 2013

### omer21

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 $u'_x(0,t)=0$.
I put
Code (Text):
T(m+1,j)=2;
just after second for loop but still i can't get correct results.

6. Mar 18, 2013

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

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

7. Mar 18, 2013

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

8. Mar 19, 2013

### jfgobin

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

I will write some code later.

9. Mar 19, 2013

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

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:

File size:
18.4 KB
Views:
180
• ###### heattrans005.png
File size:
20.3 KB
Views:
203
10. Mar 19, 2013