Matlab create finite difference matrix for Backward Euler Method

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?

Answer 1

[Torsten]

n = 98;
full(gallery('tridiag',n,-s,1+2*s,-s))

Best wishes

Torsten.

[Afzal Ali]

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.

Answer 2

[Torsten]

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.