C..2 1次元熱方程式陰解法

差分法の部分を陰解法 (いわゆる $\theta$ 法) に置き換えたプログラム heat1d_i.jl を紹介する (可視化のやり方については heat1d_e.jl と同じである)。

この項も、アルゴリズムの詳細については、「発展系の数値解析」などを参照すること。

2つのプログラムを示す。まずは LU 分解を自前の関数で実現したもの。

# heat1d_i.jl --- 1次元熱伝導方程式の初期値境界値問題を差分法(陰解法)で解く
#=  Julia v1.3 で動作確認 (2019/12/29)

 curl -O -O http://nalab.mind.meiji.ac.jp/~mk/misc/20191229/heat1d_i_v2.jl

 using Printf
 using Plots
 gr()
 include("heat1d_i.jl")
 heat1d_i(0.5, 0.01, 100, 0.5)
=#

function u0(x)
  if x < 0.5
    x
  else
    1.0-x
  end
end

# 三重対角行列のLU分解
function trilu1(al, ad, au)
  n=size(ad)[1]
  for i=1:n-1
      al[i+1] /= ad[i]
      ad[i+1] -= au[i] * al[i+1]
  end
end

# LU分解済みの三重対角行列を元に三項方程式を解く
function trisol1(al,ad,au,b)
  n=size(ad)[1]
  for i=1:n-1
    b[i+1] -= b[i] * al[i+1]
  end
  b[n] /= ad[n]
  for i=n-1:-1:1
    b[i] = (b[i] - au[i] * b[i + 1]) / ad[i]
  end
end

function heat1d_i(tMax=0.5, dt=0.01, N=100, lambda=0.5)
# tMax=0.5; dt=0.01; N=100; lambda=0.5
  h=1.0/N
  tau=lambda*h*h
  nMax = Int(ceil(tMax / tau))
  nskip = max(Int(round(dt/tau)),1)
  theta=0.5
  # 格子点
  x=range(0.0,1.0,length=N+1)
  # 初期値
  u=u0.(x)
  tstr=@sprintf("1D heat equation u_t=u_{xx}, t=%.3f",0.0)
  p=plot(x,u,title=tstr,xaxis=("x"),yaxis=("u",(-0.02,0.5)),legend=false)
  display(p)
  sleep(1)
  #
  println("N=$N, tMax=$tMax, h=$h, tau=$tau, lambda=$lambda, nMax=$nMax")
  # 連立1次方程式の係数行列
  ad = (1.0+2.0*theta*lambda)*ones(N-1)
  al = - theta*lambda*ones(N-1)
  au = - theta*lambda*ones(N-1)
  trilu1(al,ad,au)
  # 連立1次方程式の右辺の準備
  c3=1.0-2.0*(1.0-theta)*lambda
  c4=(1.0-theta)*lambda
  # ベクトルの用意 (不要かも)
  # b=zeros(N-1)
  newu=similar(u)
  for n=1:nMax
    t=n*tau
    b=c3*u[2:end-1]+c4*(u[1:end-2]+u[3:end])
    trisol1(al,ad,au,b)
    u[2:end-1]=b
    u[1] = u[end] = 0.0 # 同次Dirichlet境界条件
    # println(u)
    if (n % nskip == 0 || n == nMax)
      tstr=@sprintf("1D heat equation u_t=u_{xx}, t=%.3f",t)
      p=plot(x,u,title=tstr,xaxis=("x"),yaxis=("u",(-0.02,0.5)),legend=false)
      display(p)
      sleep(0.01)
    end
  end
end

実は、heat1d_i.jl では、最初 al = au = -theta*lambda*ones(N-1) としていて、うまく動かず、原因に気づくまで苦労した。 al と au は最初は同じ値を持つが、 LU分解した後は値が変わるはずで、 al を au の参照にしてはいけない。上のプログラムのように、個別に -theta*lambda*ones(N-1) を代入するとか、 au = copy(al) と明示的にコピーするとか、すべきである。

書き始めてしばらくは、三重対角行列のLU分解をする関数の使い方が探せなかったので、 heat1d_i.jl では、自前で用意した trulu1(), trusol1() で、その処理をしている (これらはC言語で書いたコードを移植したものである)。その後、LU 分解を標準ライブラリィに含まれる LinearAlgebra の中の Tridiagonal() を発見して作ったのが、次の heat11d_i_v2.jl である。

# heat1d_i_v2.jl --- 1次元熱伝導方程式の初期値境界値問題を差分法(陰解法)で解く
#=  Julia v1.3 で動作確認 (2020/1/2)

 curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20191229/heat1d_i_v2.jl

 using Printf
 using Plots
 gr()
 using LinearAlgebra
 include("heat1d_i_v2.jl")
 heat1d_i(0.5, 0.01, 100, 0.5)
=#

function u0(x)
  if x < 0.5
    x
  else
    1.0-x
  end
end

function heat1d_i_v2(tMax=0.5, dt=0.01, N=100, lambda=0.5)
# tMax=0.5; dt=0.01; N=100; lambda=0.5
  h=1.0/N
  tau=lambda*h*h
  nMax = Int(ceil(tMax / tau))
  nskip = max(Int(round(dt/tau)),1)
  theta=0.5
  # 格子点
  x=range(0.0,1.0,length=N+1)
  # 初期値
  u=u0.(x)
  tstr=@sprintf("1D heat equation u_t=u_{xx}, t=%.3f",0.0)
  p=plot(x,u,title=tstr,xaxis=("x"),yaxis=("u",(-0.02,0.5)),legend=false)
  display(p)
  sleep(1)
  #
  println("N=$N, tMax=$tMax, h=$h, tau=$tau, lambda=$lambda, nMax=$nMax")
  # 連立1次方程式の係数行列
  ad = (1.0+2.0*theta*lambda)*ones(N-1)
  al = - theta*lambda*ones(N-2)
  au = - theta*lambda*ones(N-2)
  a = Tridiagonal(al,ad,au)
  L,U=lu(a)
  # 連立1次方程式の右辺の準備
  c3=1.0-2.0*(1.0-theta)*lambda
  c4=(1.0-theta)*lambda
  # ベクトルの用意 (不要かも)
  #b=zeros(N-1)
  for n=1:nMax
    t=n*tau
    b=c3*u[2:end-1]+c4*(u[1:end-2]+u[3:end])
    u[2:end-1]=U\(L\b)
    u[1] = u[end] = 0.0 # 同次Dirichlet境界条件
    # println(u)
    if (n % nskip == 0 || n == nMax)
      tstr=@sprintf("1D heat equation u_t=u_{xx}, t=%.3f",t)
      p=plot(x,u,title=tstr,xaxis=("x"),yaxis=("u",(-0.02,0.5)),legend=false)
      display(p)
      sleep(0.01)
    end
  end
end

ターミナルで

$ curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20191229/heat1d_i_v2.jl $ julia julia> using Printf julia> using Plots julia> gr() julia> using LinearAlgebra julia> include("heat1d_i_v2.jl") julia> heat1d_i_v2(0.5, 0.01, 100, 0.5)

桂田祐史

C..2 1次元熱方程式 陰解法

C..2 1次元熱方程式陰解法