2次元の正方形領域 における熱方程式の 初期値境界値問題
( , ) | ||
( , ) | ||
( ) |
アルゴリズムは標準的なもので、どこでも見られるが、例えば 「熱方程式に対する差分法 I」 に説明がある。
まず可視化を gnuplot で行うバージョンを示す。
# heat2d_e_gnuplot.jl --- 空間2次元熱方程式の初期値境界値問題 (Dirichlet,Euler) #= using Printf # @printf using LinearAlgebra # norm() include("heat2d_e_gnuplot.jl") heat2d_e_gnuplot(0.05,0.0005,100,50,0.5) 80.647537 seconds (18.70 M allocations: 442.070 MiB, 0.12% gc time) =# # 初期値 function u0(x,y) sin(5 * pi * x) * sin(3 * pi * y) end # gnuplot で u を描画 function drawgraph(t, x, y, u, gp) m,n=size(u) @printf(gp, "splot '-' with lines title \"t=%.3f\"\n", t) for i=1:m for j=1:n @printf(gp, "%f %f %f\n", x[i], y[j], u[i,j]) end @printf(gp, "\n") end @printf(gp, "e\n") flush(gp) end function heat2d_e_gnuplot(tMax=0.05, dt=0.0005, Nx=100, Ny=50, lambda=0.5) a = 0.0; b = 2.0; c = 0.0; d = 1.0 hx = (b - a) / Nx; hy = (d - c) / Ny; hx2 = hx * hx; hy2 = hy * hy tau = lambda * (hx2 * hy2) / (hx2 + hy2) lambdax = tau / hx2; lambday = tau / hy2 nMax = Int(ceil(tMax / tau)) nskip = max(Int(round(dt / tau)),1) # println("tMax=$tMax, Δt=$dt, Nx=$Nx, Ny=$Ny, λ=$lambda") println("hx=$hx, hy=$hy, tau=$tau") println("減衰率=e^{-2π^2 τ}=$(exp(-2*pi*pi*tau))") println("lambda=$lambda, lambdax+lambday=$(lambdax+lambday)") println("nskip=$nskip") # メモリー確保 x = range(a, b, length = Nx+1) y = range(c, d, length = Ny+1) u = u0.(x*ones(Ny+1)', ones(Nx+1)*y') newu=zeros(size(u)) # gnuplot の起動 println("gnuplot") p=open(pipeline(`gnuplot --persist`; stderr=Pipe()), "r+") println(p.in, "set xlabel \"x\"") println(p.in, "set ylabel \"y\"") println(p.in, "set zlabel \"u\"") println(p.in, "set title \"2D heat equation u_t=u_{xx}+u_{yy}\" at screen 0.4,0.95"); println(p.in, "set xrange [a:b]") println(p.in, "set yrange [c:d]") println(p.in, "set zrange [-1.0:1.0]") println(p.in, "set hidden3d") println(p.in, "set contour") # 初期値 t=0.0 drawgraph(t,x,y,u,p.in) @printf("||u(・,%g)||=%g\n", t, norm(u)) sleep(3) # 時間発展 for n=1:nMax t=n*tau newu[2:end-1,2:end-1]=((1.0-2.0*lambda)*u[2:end-1,2:end-1] + lambdax*(u[1:end-2,2:end-1]+u[3:end,2:end-1]) + lambday*(u[2:end-1,1:end-2]+u[2:end-1,3:end])) u=copy(newu) @printf("||u(・,%g)||=%g\n", t, norm(u)) if n % nskip == 0 drawgraph(t,x,y,u,p.in) end end close(p) end
ターミナルで |
$ curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20191229/heat2d_e_gnuplot.jl $ julia julia> using Printf julia> using LinearAlgebra julia> include("heat2d_e_gnuplot.jl") julia> heat2d_e_gnuplot() |
差分方程式を
newu[2:end-1,2:end-1]=((1.0-2.0*lambda)*u[2:end-1,2:end-1] + lambdax*(u[1:end-2,2:end-1]+u[3:end,2:end-1]) + lambday*(u[2:end-1,1:end-2]+u[2:end-1,3:end])) u=copy(newu) |
for i=2:Nx for j=2:Ny newu[i,j] = ((1.0-2.0*lambda) * u[i,j] + lambdax * (u[i+1,j]+u[i-1,j]) + lambday * (u[i,j+1]+u[i,j-1])) end end for i=1:Nx+1 for j=1:Ny+1 u[i,j]=newu[i,j] end end |
次に Plots を利用したバージョン。
ターミナルで |
$ curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20200101/heat2d_e_plots.jl $ julia julia> using Printf julia> using Plots julia> gr() julia> using LinearAlgebra julia> include("heat2d_e_plots.jl") julia> heat2d_e_plots() |
plot() に渡すとき、 u を転置した u' を渡していることに注意 (' は Hermite 共役を意味するので、 transpose(u) とするのが本当かもしれない)。 転置するのが嫌ならば、初期値の設定を
u=u0.(ones(Ny+1)*x',y*ones(Nx+1)') |
newu[2:end-1,2:end-1]=((1.0-2.0*lambda)*ut[2:end-1,2:end-1] + lambdax*(ut[2:end-1,1:end-2]+ut[2:end-1,3:end]) + lambday*(ut[1:end-2,2:end-1]+ut[3:end,2:end-1])) u=copy(newu) |