Matlab create finite difference matrix for Backward Euler Method

80 views (last 30 days)
Afzal Ali
Afzal Ali on 4 Nov 2016
Edited: Torsten on 14 Nov 2016
I am trying to create a finite difference matrix to solve the 1-D heat equation (Ut = kUxx) using the backward Euler Method.
I have derived the finite difference matrix, A:
u(t+1) = inv(A)*u(t) + b, where u(t+1) u(t+1) is a vector of the spatial temperature distribution at a future time step, and u(t) is the distribution at the current time step.
The matrix A is an (n-2)-by-(n-2) matrix, where n is the size of the 1-D mesh.
u1 u2 u3 u4 . . . u_n-3 u_n-1 u_n-2
A = [1+2s -s 0 0 . . . 0 0
-s 1+2s -s 0 . . . 0 0
0 -s 1+2s -s . . . 0 0
. . . . . . . . .
. . . . . . . -s 1+2s -s
. . . . . . . 0 -s 1+2s]
I have generated this matrix using a loop, realizing that for odd rows, odd columns are 1+2s and even columns are -s, while for even rows, the opposite is true. For row i, columns <= i-2 and columns >= i + 2 are zero. The code is
L = 2; % Distance
N = 100; % Number of grid points
nu = 0.2;
mesh = linspace(0,L,N);
dx = mesh(2) - mesh(1);
tsteps = 1000; % Number of time steps
tend = 2;
t = linspace(0, tend, tsteps);
s = nu*dx^2/dt^2;
for i = 1:N-2
for j = 1:N-2
if j <= i - 2
matrix(i,j) = 0;
elseif j >= i + 2
matrix(i,j) = 0;
else
if mod(i,2) ~= 0
if mod(j,2) ~= 0
matrix(i,j) = 1+2*s;
else
matrix(i,j) = -s;
end
else
if mod(j,2) ~= 0
matrix(i,j) = -s;
else
matrix(i,j) = 1+2*s;
end
end
end
end
end
This is derived from theory presented in https://espace.library.uq.edu.au/view/UQ:239427/Lectures_Book.pdf (page 23)
Is there a shorter way to create this matrix?

Accepted Answer

Torsten
Torsten on 4 Nov 2016
n = 98;
full(gallery('tridiag',n,-s,1+2*s,-s))
Best wishes
Torsten.
  1 Comment
Afzal Ali
Afzal Ali on 12 Nov 2016
I am solving the heat equation using Backward Euler method and central difference spatial discretisation. The problem has periodic boundary conditions, so I need to use a cyclic-tridiagonal matrix:
u1 u2 u3 u4 . . . u_n-1 u_n
A = [-2k k 0 0 . . . 0 1
k -2k k 0 . . . 0 0
0 k -2k k . . . 0 0
. . . . . . . .
. . . . . . k -2k k
1 . . . . . 0 k -2k]
where k=1/dx^2
the domain is discretised with N points in space from 0 to L (length of rod):
mesh = linspace(0, L, N);
and with mesh width:
dx = mesh(2)-mesh(1);
the time steps are obtained as:
t = linspace(0, 0.1, N); % no. of time steps = no. of points in mesh
and the time step size:
dt = t(2)-t(1);
The finite difference matrix is generated as:
A = full(gallery('tridiag',N,k,-2*k,k));
A(1,end) = 1;
A(end,1) = 1;
the initial condition, ut0, is:
ut0 = sin((n*pi/L)*mesh).';
where n=2 and L=1
I then initialise an array to store the temperatures at each time step:
u_be = zeros(N,N);
and set the first column to the initial conditions:
u_be(:,1) = ut0;
I have implemented the Backward Euler as:
for i = 2:N
x = eye(N)-dt*A;
u_be(:,i) = x\u_be(:,i-1);
end
However, this does not give me the correct temperature at u0 and un, which I know from the theoretical solution (sin(sqrt(n*pi/L)*x)e^(-t*n*pi/L)) to be equal to zero. Moreover, the answer I get at u0 and un are opposite in sign. However, with periodic boundary conditions, I would expect them to be equal.
However, if I do not solve for the boundary values within the matrix, but do a separate operation like:
u_be(1,i) = (u_be(1,i-1) + k*dt*(u_be(2,i) + u_be(end-1,i)))/(1+2*k*dt);
I get the correct boundary values. This last operation is what I would expect from the calculation with the matrix as well.

Sign in to comment.

More Answers (1)

Torsten
Torsten on 14 Nov 2016
Edited: Torsten on 14 Nov 2016
You will have to introduce two ghost points:
x_(-1) = -deltax and x_(N+1) = L + deltax
The equations for these points read
u_(-1) = u_(N-1)
u_(N+1) = u_(1)
In all other points
x_(0)=0, x_(1)=deltax ,..., x_(N-1) = L-deltax, x_(N) = L
you can use the usual discretization of the heat equation.
If you include the equations for the ghost points in your system matrix (which is easiest to program), you will get a system matrix of size (N+3)x(N+3).
Best wishes
Torsten.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!