📂ジュリア

ジュリアにおける多変数関数のブロードキャス팅

概要

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