C..1 1次元熱方程式 陽解法

アルゴリズムは良く知られたものを使っている (例えば 「発展系の数値解析」)。

慣れている人向けに、差分方程式だけ書いておく (これを見ればプログラムを納得できるかも)。

$$U_{i}^{n+1} =(1-2\lambda)U_{i}^n+\lambda\left(U_{i-1}^n+U_{i+1}^n\right) i=1,2,\cdots,N-1 ,$$

$$U_{0}^{n+1}=U_{N}^{n+1}=0$$

グラフを描くために、 E 節で説明する Plots を用いている。

# heat1d_e.jl --- 1次元熱伝導方程式の初期値境界値問題を差分法(陽解法)で解く
#  Julia v1.3 で動作確認 (2019/12/29)
# julia> using Printf
# julia> using Plots
# julia> gr()
# julia> include("heat1d_e.jl")
# julia> heat1d_e(0.5, 0.01, 100, 0.5)

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

function heat1d_e(tMax,dt,N,lambda)
# u0 は関数
  h=1.0/N
  tau=lambda*h*h
  nMax = Int(ceil(tMax / tau))
  nskip = max(Int(round(dt/tau)),1)
  # 格子点
  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(3)
  #
  println("N=$N, tMax=$tMax, h=$h, tau=$tau, lambda=$lambda, nMax=$nMax")
  #
  newu=similar(u)
  for n=1:nMax
    t=n*tau
    newu[2:end-1]=(1-2*lambda)*u[2:end-1]+lambda*(u[1:end-2]+u[3:end])
    newu[1] = newu[end] = 0.0
    u=newu
    # 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

ループの中では、単に plot() を実行しても、グラフは表示されない。 明示的に display() する必要がある。

時刻ごとのグラフを比較できるように、 座標軸の範囲が変わらないようにすることが必要で、 yaxis=(...,(-0.02,0.5),...) としたのは大事なところである。

ターミナルで

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