アルゴリズムは良く知られたものを使っている (例えば 「発展系の数値解析」)。
慣れている人向けに、差分方程式だけ書いておく (これを見ればプログラムを納得できるかも)。
( ) | ||
グラフを描くために、 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) |