ジュリアにおける多変数関数のブロードキャス팅
概要
Juliaで多変数関数をブロードキャストする方法を紹介する。Pythonなどで行うように、meshgridを作成する方法もあるし、各次元ごとにベクトルを作成して簡単に計算することもできる。
2変数関数
$$ u(t,x) = \sin(\pi x) e^{-\pi^{2}t} $$
上のような関数を$(t,x) \in [0, 0.35] \times [-1,1]$でプロットしたい場合、次のように関数値を計算できる。
x = LinRange(-1., 1, 100)
t = LinRange(0., 0.35, 200)'
u1 = @. sin(π*x)*exp(- π^2 * t)
heatmap(t', x, u1, xlabel="t", ylabel="x", title="Fig. 1")
関数自体を定義し、次のように2次元グリッドを作成して同じ結果を得ることができる。
U(t,x) = sin(π*x)*exp(- π^2 * t)
x = LinRange(-1., 1, 100)
t = LinRange(0., 0.35, 200)'
X = x * fill!(similar(t), 1)
T = fill!(similar(x), 1) * t
u2 = U.(T,X)
heatmap(t', x, u2, xlabel="t", ylabel="x", title="Fig. 2")
3変数関数
$$ u(x,y,t) = e^{-x^{2} - 2y^{2}}e^{-\pi^{2}t} $$
時空間ドメイン$(x,y,t) \in [-1,1] \times [-1,1] \times [0, 0.35]$で$u$の関数値を得たい場合、各変数の次元にのみサイズがあるようにベクトルを作成してブロードキャストすればいい。
3次元メッシュを作成してブロードキャストしたい場合は、ここを参照。
julia> x = reshape(LinRange(-1., 1, 100), (100,1,1))
100×1×1 reshape(::LinRange{Float64, Int64}, 100, 1, 1) with eltype Float64:
julia> y = reshape(LinRange(-1., 1, 100), (1,100,1))
1×100×1 reshape(::LinRange{Float64, Int64}, 1, 100, 1) with eltype Float64:
julia> t = reshape(LinRange(0.,0.35, 200), (1,1,200))
1×1×200 reshape(::LinRange{Float64, Int64}, 1, 1, 200) with eltype Float64:
julia> u3 = @. exp(-x^2) * exp(-2y^2) * exp(- π^2 * t)
100×100×200 Array{Float64, 3}:
anim = @animate for i ∈ 1:200
surface(u3[:,:,i], zlims=(0,1), clim=(-1,1))
end
コード詳細
using Plots
cd = @__DIR__
# Fig. 1
x = LinRange(-1., 1, 100)
t = LinRange(0., 0.35, 200)'
u1 = @. sin(π*x)*exp(- π^2 * t)
heatmap(t', x, u1, xlabel="t", ylabel="x", title="Fig. 1")
savefig(cd*"/fig1.png")
# Fig. 2
U(t,x) = sin(π*x)*exp(- π^2 * t)
x = LinRange(-1., 1, 100)
t = LinRange(0., 0.35, 200)'
X = x * fill!(similar(t), 1)
T = fill!(similar(x), 1) * t
u2 = U.(T,X)
heatmap(t', x, u2, xlabel="t", ylabel="x", title="Fig. 2")
savefig(cd*"/fig2.png")
# gif 1
x = reshape(LinRange(-1., 1, 100), (100,1,1))
y = reshape(LinRange(-1., 1, 100), (1,100,1))
t = reshape(LinRange(0.,0.35, 200), (1,1,200))
u3 = @. exp(-x^2) * exp(-2y^2) * exp(- π^2 * t)
anim = @animate for i ∈ 1:200
surface(u3[:,:,i], zlims=(0,1), clim=(-1,1), title="Anim. 1")
end
gif(anim, cd*"/anim1.gif", fps=30)
環境
- OS: Windows11
- Version: Julia v1.8.3, Plots v1.38.6