줄리아에서 다변수함수의 브로드캐스팅
개요
줄리아에서 다변수함수를 브로드캐스팅하는 방법을 소개한다. 파이썬 등에서 하는 것처럼 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