# 総和がNのN個の正の実数を成分とするランダムベクトルの生成
using Random: AbstractRNG, SamplerTrivial, rand!, default_rng
struct Simplex{T} N::Int end
Simplex(N) = Simplex{Float64}(N)
function Base.rand(rng::AbstractRNG, d::SamplerTrivial{Simplex{T}}) where T
N = d[].N
X = Vector{T}(undef, N)
rand!(rng, @view(X[begin:end-1]))
X[end] = 1
X .*= N
sort!(X)
@. @views X[begin+1:end] -= X[begin:end-1]
X
end
rand(Simplex(3))
3-element Vector{Float64}: 0.3821729604174189 1.6696744406662767 0.9481525989163044
rand(Simplex(10)) |> sum
10.0
rand(Simplex(3), 5)
5-element Vector{Any}: [0.8831689667896425, 2.0525430221737353, 0.06428801103662218] [0.8351431822542357, 1.4112466595378925, 0.7536101582078718] [0.9081146178891397, 2.0247026093705247, 0.06718277274033557] [0.2521109775292505, 0.9398980255839358, 1.8079909968868138] [1.0304935927256627, 1.6274943685203813, 0.342012038753956]
# 総和がNのN個の正の実数を成分とするランダムベクトルの成分の分布は
# Nを大きくすると期待値1の指数分布に近付く.
# これは「統計力学」の最も易しい場合になっている.
using Plots
N = 10^4
histogram(rand(Simplex(N)); norm=true, alpha=0.3, label="rand(Simplex($N))")
plot!(x->exp(-x), 0, 8; label="exp(-x)", lw=2)
plot!(xlim=(0, 8))
N = length(X)が大きなとき, Xのすべての成分を正に保ったままで, N == sum(X)が不変になるようなランダムウォークを繰り返すと, Xの成分の分布は期待値1の指数分布で近似されるようになる.
# 富のランダム再分配
function randomly_redistribute!(rng::AbstractRNG, X, h)
N = length(X)
i = rand(rng, 1:N)
j = rand(rng, 1:N-1)
j = ifelse(j == i, N, j)
dx = h*abs(randn())
if dx < X[i]
X[i] -= dx
X[j] += dx
end
X
end
function iter_redistribute!(X, L, h)
rng = default_rng()
for _ in 1:L
randomly_redistribute!(rng, X, h)
end
X
end
iter_redistribute! (generic function with 1 method)
N = 10^4
X = ones(N)
L = 10^6
h = 0.02
tmax = 200
anim = @animate for t in [fill(0, 10); 1:tmax-1; fill(tmax, 20)]
t == 0 || t == tmax || iter_redistribute!(X, L, h)
histogram(X; norm=true, alpha=0.3, label="MCMC")
plot!(x -> exp(-x), 0, 8; label="exp(-x)", lw=2)
plot!(xlim=(-0.1, 8), ylim=(-0.02, 1.05))
title!("t = $t"; titlefontsize=12)
end
gif(anim, "exp1.gif", fps=20)
┌ Info: Saved animation to │ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\exp1.gif └ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
N = length(X)が大きなとき, Xのすべての成分が正でかつ, sum(X)がN以下であるという条件を保ったままランダムウォークを繰り返すと, Xの成分の分布は期待値1の指数分布で近似されるようになる. sum(X)をぴったりNに固定せずに, N以下に広げても, 収束する先の分布は変わらない.
詳しい一般論については, 『Kullback-Leibler情報量とSanovの定理』の第4.2節を参照せよ. そこではエネルギー $E_i$ の平均値をぴったり固定せずに不等号で制限している.
function iter_randwalk!(X, L, h)
N = length(X)
S = sum(X)
if S > N
X ./= 1.1S
end
rng = default_rng()
for _ in 1:L
i = rand(rng, 1:N)
dx = h*randn()
if X[i] + dx > 0 && S + dx ≤ N
X[i] += dx
S += dx
end
end
X
end
iter_randwalk! (generic function with 1 method)
N = 10^4
X = 0.8ones(N) # テキトーな初期条件
L = 10^6
h = 0.02
tmax = 200
anim = @animate for t in [fill(0, 10); 1:tmax-1; fill(tmax, 20)]
t == 0 || t == tmax || iter_randwalk!(X, L, h)
histogram(X; norm=true, alpha=0.3, label="MCMC")
plot!(x -> exp(-x), 0, 8; label="exp(-x)", lw=2)
plot!(xlim=(-0.1, 8), ylim=(-0.02, 1.05))
title!("t = $t"; titlefontsize=12)
end
gif(anim, "exp2.gif", fps=20)
┌ Info: Saved animation to │ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\exp2.gif └ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
# Euclidノルムの二乗がNのランダムベクトルを生成
using Random: AbstractRNG, SamplerTrivial, rand!, default_rng
using LinearAlgebra: norm, dot
norm2(x) = dot(x, x)
struct Sphere{T} N::Int end
Sphere(N) = Sphere{Float64}(N)
function Base.rand(rng::AbstractRNG, d::SamplerTrivial{Sphere{T}}) where T
N = d[].N
X = randn(rng, T, N)
X .*= √N/norm(X)
X
end
X = rand(Sphere(4))
4-element Vector{Float64}: -1.784652385724709 0.8551989101001829 -0.06385898611929577 -0.2820863629854841
X = rand(Sphere(10)) |> norm2
10.000000000000002
rand(Sphere(3), 5)
5-element Vector{Any}: [-0.9650572142034735, 1.1297418111946969, 0.8901393224391091] [-1.043613527473952, -0.9064699885083874, 1.0436392888383315] [0.3038286084030116, 0.018722254709423625, 1.7050916848939488] [1.2764177438394357, 0.9505113380289749, -0.6835831620879875] [0.6999735988769156, 1.3721366361593437, -0.7920088462792703]
# Euclidノルムの二乗がNのランダムベクトルの成分の分布は
# Nを大きくすると標準正規分布に近付く.
# これは「統計力学」におけるMaxwell–Boltzmann分布の場合になっている.
using Plots
N = 10^4
histogram(rand(Sphere(N)); norm=true, alpha=0.3, label="rand(Sphere($N))")
plot!(x->exp(-x^2/2)/√(2π), -4, 4; label="standard normal dist.", lw=2)
N = length(X)が大きなとき, XのEuclidノルムの2乗がN以下であるという条件を保ったままランダムウォークを繰り返すと, Xの成分の分布は標準正規分布で近似されるようになる. XのEuclidノルムの2乗をぴったりNに固定せずに, N以下に広げても, 収束する先の分布は変わらない.
詳しい一般論については, 『Kullback-Leibler情報量とSanovの定理』の第4.2節を参照せよ. そこではエネルギー $E_i$ の平均値をぴったり固定せずに不等号で制限している.
function iter_randwalk_in_ball!(X, L, h)
N = length(X)
E = norm2(X) # E stands for the total Energy
if E > N
X ./= 1.1*√E
end
rng = default_rng()
for _ in 1:L
i = rand(rng, 1:N)
dx = h*randn()
dE = 2*X[i]*dx + dx^2
if E + dE ≤ N
X[i] += dx
E += dE
end
end
X
end
iter_randwalk_in_ball! (generic function with 1 method)
N = 10^4
X = [0.8 .+ 0.5rand(N ÷ 2); -0.9 .+ 0.2randn(N ÷ 2)] # テキトーな初期条件
L = 10^6
h = 0.01
tmax = 200
anim = @animate for t in [fill(0, 10); 1:tmax-1; fill(tmax, 20)]
t == 0 || t == tmax || iter_randwalk_in_ball!(X, L, h)
histogram(X; norm=true, alpha=0.3, label="MCMC")
plot!(x->exp(-x^2/2)/√(2π), -4, 4; label="standard normal dist.", lw=2)
plot!(xlim=(-4, 4), ylim=(-0.02, 0.5))
title!("t = $t"; titlefontsize=12)
end
gif(anim, "normal2.gif", fps=20)
┌ Info: Saved animation to │ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\normal2.gif └ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
N = length(X)が大きなとき, Xの各成分が正で加法平均がμ以下であるという条件だけではなく, 対数平均がある値ν以上になるという条件を保ったままランダムウォークを繰り返すと, Xの成分の分布はガンマ分布で近似されるようになる.
using Random: default_rng
using Plots
using SpecialFunctions
using Distributions
using Printf
rd(x, d=2) = round(x; digits=d)
using ProgressMeter
function iter_randwalk_gamma!(X, μ, ν, L, h)
N = length(X)
S = mean(X)
U = mean(log, X)
rng = default_rng()
for _ in 1:L
i = rand(rng, 1:N)
dx = h*randn()
X_i_new = X[i] + dx
if X_i_new > 0
S_new = S + dx/N
U_new = U + (-log(X[i]) + log(X_i_new))/N
if S_new ≤ μ && U_new ≥ ν
X[i] = X_i_new
S = S_new
U = U_new
end
end
end
X
end
Delta(α) = log(α) - digamma(α)
Delta (generic function with 1 method)
plot(Delta, 0.2, 10; label="Delta(α)", size=(400, 250))
function anim_gamma!(μ, α, X, L, h, tmax, xlim=(0, 10), ylim=(0, 0.6))
θ = μ/α
Δ = Delta(α)
ν = log(μ) - Δ
S = [mean(X)]
U = [mean(log, X)]
prog = Progress(tmax+30, 0)
anim = @animate for t in [fill(0, 10); 1:tmax; fill(tmax+1, 20)]
if 1 ≤ t ≤ tmax
iter_randwalk_gamma!(X, μ, ν, L, h)
push!(S, mean(X))
push!(U, mean(log, X))
end
P = histogram(X; norm=true, alpha=0.3, label="MCMC")
title!("μ = $(rd(μ)), ν = log(μ) - $(rd(Δ,4)), t = $t")
f_label = @sprintf "Gamma(%.2f, %.2f)" α θ
plot!(x -> pdf(Gamma(α, θ), x), xlim...; label=f_label, color=:red, lw=2)
plot!(; xlim, ylim)
Q = hline([μ]; label="", color=:red, lw=2)
plot!(0:min(tmax,t), S; label="", color=:blue)
title!("mean(X)")
plot!(xlim=(0, tmax), ylim=(minimum(S), μ + 0.05*(μ - minimum(S))))
R = hline([ν]; label="", color=:red, lw=2)
plot!(0:min(tmax,t), U; label="", color=:blue)
title!("mean(log, X)")
plot!(xlim=(0, tmax), ylim=(ν - 0.05*(maximum(U) - ν), maximum(U)))
layout = @layout[a [b; c]]
plot(P, Q, R; layout, size=(800, 300))
plot!(legendfontsize=8, titlefontsize=10)
next!(prog)
end
gif(anim, @sprintf("Gamma(%.2f, %.2f).gif", α, θ), fps=20)
end
anim_gamma! (generic function with 3 methods)
μ = 2.0
α = 1.0
N = 10^4
X = 0.7μ*ones(N)
L = 10^6
h = 0.05
tmax = 200
xlim=(0, 10)
ylim=(0, 0.6)
anim_gamma!(μ, α, X, L, h, tmax, xlim, ylim)
Progress: 100%%|█████████████████████████████████████████| Time: 0:00:28
┌ Info: Saved animation to
│ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\Gamma(1.00, 2.00).gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
μ = 2.0
α = 5.0
N = 10^4
X = 0.95μ*ones(N)
L = 8*10^5
h = 0.01
tmax = 200
xlim = (0, 8)
ylim = (0, 0.6)
anim_gamma!(μ, α, X, L, h, tmax, xlim, ylim)
Progress: 100%%|█████████████████████████████████████████| Time: 0:00:23
┌ Info: Saved animation to
│ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\Gamma(5.00, 0.40).gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
μ = 2.0
α = 25.0
N = 10^4
X = 0.99μ*ones(N)
L = 5*10^5
h = 0.005
tmax = 200
xlim = (0, 5)
ylim = (0, 1.2)
anim_gamma!(μ, α, X, L, h, tmax, xlim, ylim)
Progress: 100%%|█████████████████████████████████████████| Time: 0:00:20
┌ Info: Saved animation to
│ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\Gamma(25.00, 0.08).gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
using QuadGK
using NLsolve
bf(x, a, f=√) = exp(-a[1]*x + (a[2]-1)*f(x))
Zf(a, f=√, xmax=1e3) = quadgk(x -> bf(x, a, f), 0, xmax)[1]
function Eq!(p, a, f=√, μ=2.0, ν=√μ - 0.05, xmax=1e3)
Z = Zf(a, f, xmax)
p[1] = quadgk(x -> x*bf(x, a, f), 0, xmax)[1]/Z - μ
p[2] = quadgk(x -> f(x)*bf(x, a, f), 0, xmax)[1]/Z - ν
p
end
function solve_Eq(f=√, μ=2.0, ν=√μ - 0.05, xmax=1e3)
F!(p, a) = Eq!(p, a, f, μ, ν, xmax)
sol = nlsolve(F!, ones(2))
a = sol.zero
Z = Zf(a, f, xmax)
p(x) = bf(x, a, f)/Z
(a = a, p = p)
end
solve_Eq (generic function with 5 methods)
@time sol = nlsolve(Eq!, ones(2))
4.582457 seconds (4.38 M allocations: 244.401 MiB, 2.42% gc time)
Results of Nonlinear Solver Algorithm * Algorithm: Trust-region with dogleg and autoscaling * Starting Point: [1.0, 1.0] * Zero: [3.2418318581941707, 9.039303922174012] * Inf-norm of residuals: 0.000000 * Iterations: 8 * Convergence: true * |x - x'| < 0.0e+00: false * |f(x)| < 1.0e-08: true * Function Calls (f): 9 * Jacobian Calls (df/dx): 9
↑コンパイルによる遅延時間が含まれている.
数値積分を3回実行して計算される函数の零点を見付ける問題を解かせることになるので, 実行時間が長くなり過ぎることを心配したが, Iterations: 8 なので大丈夫そうである.
dump(sol)
NLsolve.SolverResults{Float64, Float64, Vector{Float64}, Vector{Float64}} method: String "Trust-region with dogleg and autoscaling" initial_x: Array{Float64}((2,)) [1.0, 1.0] zero: Array{Float64}((2,)) [3.2418318581941707, 9.039303922174012] residual_norm: Float64 7.472960028565012e-10 iterations: Int64 8 x_converged: Bool false xtol: Float64 0.0 f_converged: Bool true ftol: Float64 1.0e-8 trace: NLsolve.SolverTrace states: Array{NLsolve.SolverState}((0,)) f_calls: Int64 9 g_calls: Int64 9
@time solve_Eq()
0.256688 seconds (714.94 k allocations: 40.283 MiB, 6.13% gc time)
(a = [3.2418318581941707, 9.039303922174012], p = var"#p#20"{typeof(sqrt), Float64, Vector{Float64}}(sqrt, 356.6177646324729, [3.2418318581941707, 9.039303922174012]))
↑二度目の実行はそこそこ速い.
以下, $N$ は十分に大きいと仮定する.
$f(x)$ は正の実数 $x$ について定義された, 上に凸な狭義単調増加函数であるとする(所謂効用函数). パラメーター $a=(a_1,a_2)$ を持つ $x$ の確率密度函数 $p(x|a)$ を
$$ p(x|a) = Z(a,b)^{-1}e^{-a_1 x + (a_2 - 1)f(x)}, \quad Z(a,b) = \int_0^\infty e^{-a_1 x + (a_2 - 1)f(x)}\,dx $$と定める. $\mu$ は正の実数であるとし, $\nu$ は
$$ \frac{1}{\mu}\int_0^\infty f(x)e^{-x/\mu}\,dx = \int_0^\infty f(x)p(x|1/\mu, 1)\,dx \le \nu < \log\mu $$を満たしていると仮定する. $a=(a_1,a_2)$ をこのような $\mu, \nu$ から
$$ \int_0^\infty x\,p(x|a)\,dx = \mu, \quad \int_0^\infty f(x)p(x|a)\,dx = \nu $$という条件によって定める.
実数達 $X_i$ ($i=1,2,\ldots,N$) の範囲に以下の制限を課す:
$$ X_i > 0, \quad \frac{1}{N}\sum_{i=1}^N X_i \le \mu, \quad \frac{1}{N}\sum_{i=1}^N f(X_i) \ge \nu. $$この条件のもとで, $X=(X_1,\ldots,X_N)$ をランダムウォークさせると, $X_i$ 達のなす分布(大雑把には $X_i$ 達のヒストグラムだとみなしてよい)は, 上のように $\mu, \nu$ から $a=(a_1,a_2)$ を定めたときの確率密度函数 $p(x|a)$ を持つ分布で近似される.
このことを $f(x) = \sqrt{x}$ の場合について数値的に確認しよう.
using Random: default_rng
using Plots
using Printf
rd(x, d=2) = round(x; digits=d)
function iter_randwalk_lbfun!(f, X, μ, ν, L, h)
N = length(X)
S = mean(X)
U = mean(f, X)
rng = default_rng()
for _ in 1:L
i = rand(rng, 1:N)
dx = h*randn()
X_i_new = X[i] + dx
if X_i_new > 0
S_new = S + dx/N
U_new = U + (-f(X[i]) + f(X_i_new))/N
if S_new ≤ μ && U_new ≥ ν
X[i] = X_i_new
S = S_new
U = U_new
end
end
end
X
end
iter_randwalk_lbfun! (generic function with 1 method)
μ = 2.0
ν = √μ - 1e8
νmin = quadgk(x -> √x*bf(x, [1/μ, 1], √), 0, 1e3)[1]/Zf([1/μ, 1], √)
@show νmin
@show Δmax = √μ - νmin
N = 10^4
X = 0.99μ*ones(N)
L = 10^6
h = 0.05
tmax = 100
@time iter_randwalk_lbfun!(√, X, μ, ν, tmax*L, h)
@show mean(√, X)
@show √μ - mean(√, X)
histogram(X; norm=true, alpha=0.3, label="")
plot!(x -> exp(-x/μ)/μ, 0, 15; label="", lw=2)
νmin = 1.2533141378760635 Δmax = √μ - νmin = 0.16089942449703165 3.541441 seconds (143.91 k allocations: 8.661 MiB) mean(√, X) = 1.2581386302589952 √μ - mean(√, X) = 0.15607493211409995
μ = 2.0
ν = √μ - Δmax
N = 10^4
X = 0.99μ*ones(N)
L = 10^6
h = 0.05
tmax = 100
@time iter_randwalk_lbfun!(√, X, μ, ν, tmax*L, h)
@show mean(√, X)
@show Δmax = √μ - mean(√, X)
histogram(X; norm=true, alpha=0.3, label="")
plot!(x -> exp(-x/μ)/μ, 0, 15; label="", lw=2)
3.686462 seconds (1 allocation: 16 bytes) mean(√, X) = 1.2602651969994991 Δmax = √μ - mean(√, X) = 0.15394836537359602
μ = 2.0
ν = √μ - 0.1
@time a, p = solve_Eq(√, μ, ν)
N = 10^4
X = 0.99μ*ones(N)
L = 10^6
h = 0.05
tmax = 100
@time iter_randwalk_lbfun!(√, X, μ, ν, tmax*L, h)
histogram(X; norm=true, alpha=0.3, label="")
plot!(p, 0, 12; label="", lw=2)
0.015733 seconds (44.19 k allocations: 1.164 MiB) 3.626243 seconds (1 allocation: 16 bytes)
μ = 2.0
ν = √μ - 0.05
@time a, p = solve_Eq(√, μ, ν)
N = 10^4
X = 0.99μ*ones(N)
L = 10^6
h = 0.05
tmax = 100
@time iter_randwalk_lbfun!(√, X, μ, ν, tmax*L, h)
histogram(X; norm=true, alpha=0.3, label="")
plot!(p, 0, 8; label="", lw=2)
0.003626 seconds (49.61 k allocations: 1.152 MiB) 3.801518 seconds (1 allocation: 16 bytes)
μ = 2.0
ν = √μ - 0.02
@time a, p = solve_Eq(√, μ, ν)
N = 10^4
X = 0.99μ*ones(N)
L = 10^6
h = 0.05
tmax = 100
@time iter_randwalk_lbfun!(√, X, μ, ν, tmax*L, h)
histogram(X; norm=true, alpha=0.3, label="")
plot!(p, 0, 5; label="", lw=2)
0.003115 seconds (40.43 k allocations: 949.234 KiB) 3.814228 seconds (1 allocation: 16 bytes)
μ = 2.0
ν = √μ - 0.01
@time a, p = solve_Eq(√, μ, ν)
N = 10^4
X = 0.99μ*ones(N)
L = 10^6
h = 0.05
tmax = 100
@time iter_randwalk_lbfun!(√, X, μ, ν, tmax*L, h)
histogram(X; norm=true, alpha=0.3, label="")
plot!(p, 0, 5; label="", lw=2)
0.003681 seconds (47.33 k allocations: 1.071 MiB) 4.226219 seconds (1 allocation: 16 bytes)
以上を見ると理論通りの分布に収束していることが分かる.
using ProgressMeter
function anim_lbfun!(f, μ, Δ, X, L, h, tmax, xlim=(0, 10), ylim=(0, 0.6))
ν = f(μ) - Δ
a, p = solve_Eq(f, μ, ν)
M = [mean(X)]
U = [mean(f, X)]
prog = Progress(tmax+30, 0)
anim = @animate for t in [fill(0, 10); 1:tmax; fill(tmax+1, 20)]
if 1 ≤ t ≤ tmax
iter_randwalk_lbfun!(f, X, μ, ν, L, h)
push!(M, mean(X))
push!(U, mean(f, X))
end
P = histogram(X; norm=true, alpha=0.3, label="MCMC")
title!("μ = $(rd(μ)), ν = $f(μ) - $(rd(Δ,4)), t = $t")
f_label = @sprintf "Dist(%s, %.2f, %.2f)" f a[1] a[2]
plot!(p, xlim...; label=f_label, color=:red, lw=2)
plot!(; xlim, ylim)
Q = hline([μ]; label="", color=:red, lw=2)
plot!(0:min(tmax,t), M; label="", color=:blue)
title!("mean(X)")
plot!(xlim=(0, tmax), ylim=(minimum(M), μ + 0.05*(μ - minimum(M))))
R = hline([ν]; label="", color=:red, lw=2)
plot!(0:min(tmax,t), U; label="", color=:blue)
title!("mean($f, X)")
plot!(xlim=(0, tmax), ylim=(ν - 0.05*(maximum(U) - ν), maximum(U)))
layout = @layout[a [b; c]]
plot(P, Q, R; layout, size=(800, 300))
plot!(legendfontsize=8, titlefontsize=10)
next!(prog)
end
gif(anim, @sprintf("Dist(%s, %.2f, %.2f).gif", f, a[1], a[2]), fps=20)
end
anim_lbfun! (generic function with 3 methods)
μ = 2.0
Δ = 0.16
N = 10^4
X = 0.99μ*ones(N)
L = 5*10^5
h = 0.05
tmax = 200
xlim = (0, 15)
ylim = (0, 0.6)
anim_lbfun!(√, μ, Δ, X, L, h, tmax, xlim, ylim)
Progress: 100%%|█████████████████████████████████████████| Time: 0:00:19
┌ Info: Saved animation to
│ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\Dist(sqrt, 0.51, 1.02).gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
μ = 2.0
Δ = 0.1
N = 10^4
X = 0.99μ*ones(N)
L = 10^5
h = 0.05
tmax = 200
xlim = (0, 12)
ylim = (0, 0.5)
anim_lbfun!(√, μ, Δ, X, L, h, tmax, xlim, ylim)
Progress: 100%%|█████████████████████████████████████████| Time: 0:00:16
┌ Info: Saved animation to
│ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\Dist(sqrt, 1.31, 3.48).gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
μ = 2.0
Δ = 0.05
N = 10^4
X = 0.99μ*ones(N)
L = 10^5
h = 0.05
tmax = 200
xlim = (0, 8)
ylim = (0, 0.5)
anim_lbfun!(√, μ, Δ, X, L, h, tmax, xlim, ylim)
Progress: 100%%|█████████████████████████████████████████| Time: 0:00:16
┌ Info: Saved animation to
│ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\Dist(sqrt, 3.24, 9.04).gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
μ = 2.0
Δ = 0.02
N = 10^4
X = 0.99μ*ones(N)
L = 2*10^4
h = 0.05
tmax = 200
xlim = (0, 5)
ylim = (0, 0.7)
anim_lbfun!(√, μ, Δ, X, L, h, tmax, xlim, ylim)
Progress: 100%%|█████████████████████████████████████████| Time: 0:00:17
┌ Info: Saved animation to
│ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\Dist(sqrt, 8.62, 24.29).gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104
μ = 2.0
Δ = 0.01
N = 10^4
X = 0.99μ*ones(N)
L = 10^4
h = 0.05
tmax = 200
xlim = (0, 5)
ylim = (0, 0.9)
anim_lbfun!(√, μ, Δ, X, L, h, tmax, xlim, ylim)
Progress: 100%%|█████████████████████████████████████████| Time: 0:00:16
┌ Info: Saved animation to
│ fn = C:\Users\genkuroki\OneDrive\Math\Math0045\Dist(sqrt, 17.48, 49.36).gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\fewot\src\animation.jl:104