Matlab create finite difference matrix for Backward Euler Method
13 views (last 30 days)
Show older comments
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?
0 Comments
Accepted Answer
More Answers (1)
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.
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!