1次元の熱伝導方程式を有限要素法で解きます。
定式化
はじめに，定式化を行います。熱伝導方程式は次式で表せます。
u^\prime(x, t) = \kappa \frac{\partial^2 u}{\partial x^2}
いま，境界における温度勾配がゼロであるようなモデルを考え（自然境界条件），試験関数$\phi$を用いてこの方程式を弱形式化すると，
\int u^\prime(x, t) \phi(x) dx
=-\kappa \int \frac{\partial u}{\partial x} \frac{\partial \phi}{\partial x} dx
が得られます。さらに，時間微分に1次の前進差分を適用します（上付き添字によって時刻を表すことにする）。
\int u^{n+1}(x) \phi(x) dx = \int u^{n}(x) \phi(x) dx -\kappa \Delta t \int \frac{\partial u}{\partial x} \frac{\partial \phi}{\partial x} dx
この方程式に対して，その近似解が区分的線形関数を用いて，
u(x) = \sum_{j=0}^{N} u_j \phi_j(x)
のように表現できるとすると，以下の計算式
\int_{x_i}^{x_{i+1}} \phi_i(x) \phi_{i+1}(x) dx = \frac{h}{6},
\quad
\int_{x_{i-1}}^{x_{i+1}} \phi_i(x) \phi_i(x) dx = \frac{2h}{3},
\int_{x_i}^{x_{i+1}} \frac{\partial \phi_i}{\partial x} \frac{\partial \phi_{i+1}}{\partial x} dx = -\frac{1}{h},
\quad
\int_{x_{i-1}}^{x_{i+1}} \frac{\partial \phi_i}{\partial x} \frac{\partial \phi_i}{\partial x} dx = \frac{2}{h}
が成り立つことから（ここで$h$は区間幅を表します），次式を得ます。
\displaylines{
\frac{h}{6} u_{i-1}^{n+1} + \frac{2h}{3} u_i^{n+1} + \frac{h}{6} u_{i+1}^{n+1}
= \\ \left(\frac{h}{6} + \frac{\kappa \Delta t}{h}\right) u_{i-1}^n + \left(\frac{2h}{3} - \frac{2 \kappa \Delta t}{h}\right)u_i^n + \left(\frac{h}{6} + \frac{\kappa \Delta t}{h}\right)u_{i+1}^n
}
プログラミング
先に導いた方程式に従って左辺の係数行列と右辺のベクトルを計算し，それぞれ amat, bvec に保持します。
subroutine assemble_matrix(n, h, amat)
integer, intent(in) :: n
real(real64), intent(in) :: h
real(real64), intent(out) :: amat(n,3)
real(real64), parameter :: one6th = 1.0_real64 / 6.0_real64
real(real64), parameter :: two3rd = 2.0_real64 / 3.0_real64
integer :: i
amat(:,:) = 0.0_real64
do i = 2, n-1
amat(i, 1) = one6th * h
amat(i, 2) = two3rd * h
amat(i, 3) = one6th * h
end do
amat(1, 2) = 1.0_real64
amat(1, 3) = -1.0_real64
amat(n, 1) = -1.0_real64
amat(n, 2) = 1.0_real64
end subroutine
subroutine assemble_rhs(n, h, kappa, dt, u, bvec)
integer, intent(in) :: n
real(real64), intent(in) :: h
real(real64), intent(in) :: kappa, dt
real(real64), intent(in) :: u(n)
real(real64), intent(out) :: bvec(n)
real(real64), parameter :: one6th = 1.0_real64 / 6.0_real64
real(real64), parameter :: two3rd = 2.0_real64 / 3.0_real64
integer :: i
do i = 2, n-1
bvec(i) = (one6th * h + kappa * dt / h) * u(i-1) &
+ (two3rd * h - 2.0_real64 * kappa * dt / h) * u(i) &
+ (one6th * h + kappa * dt / h) * u(i+1)
end do
bvec(1) = 0.0_real64
bvec(n) = 0.0_real64
end subroutine
LAPACK の帯行列ソルバ "GBTRF/GBTRS" を使用して LU 分解で解きます。
subroutine solve_lu_gb(n, amat, bvec)
integer, intent(in) :: n
real(real64), intent(in) :: amat(n,3)
real(real64), intent(inout) :: bvec(n)
character(len=1), parameter :: trans = 'N'
integer, parameter :: kl = 1
integer, parameter :: ku = 1
integer, parameter :: ldab = 2 * kl + ku + 1
integer :: i, j, k, jj, nrhs, ldb, info
integer :: ipiv(n)
real(real64) :: ab(ldab,n)
nrhs = 1
ldb = n
k = kl + ku + 1
do i = 1, n
do j = max(i-kl,1), min(i+ku,n)
jj = j - i + 2
ab(k+i-j,j) = amat(i,jj)
end do
end do
call dgbtrf(n, n, kl, ku, ab, ldab, ipiv, info)
if (info /= 0) error stop ': error detected!'
call dgbtrs(trans, n, kl, ku, nrhs, ab, ldab, ipiv, bvec, ldb, info)
end subroutine
きちんと解けているようです。
