Matlab create finite difference matrix for Backward Euler Method
33 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
更多回答(1 个)
Torsten
2016-11-14
编辑:Torsten
2016-11-14
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!