Matlab create finite difference matrix for Backward Euler Method
10 ビュー (過去 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!