Matlab create finite difference matrix for Backward Euler Method
10 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Afzal Ali
am 4 Nov. 2016
Bearbeitet: Torsten
am 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?
0 Kommentare
Akzeptierte Antwort
Weitere Antworten (1)
Torsten
am 14 Nov. 2016
Bearbeitet: Torsten
am 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 Kommentare
Siehe auch
Kategorien
Mehr zu Geometry and Mesh finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!