偏微分方程式: 多次元問題(2)

空間次元が多次元の偏微分方程式(続き)

今度は、逃げる・追いかける効果をモデルに入れてみよう

前回はあくまで計算に慣れるという目的のもと、おとなしめなモデルを扱った．今回は、被食者・捕食者問題の性質をよりよく反映させるため、「逃げる・追いかける」効果を先のモデルに追加しよう．

これは丁寧に考えるとそう難しくなく(授業時に説明する)、

魚 $u(x,t)$ がサメ $v(x,t)$ から逃げる効果: $du/dt$ に $\nabla ( u \nabla v)$ を足す
サメ $v(x,t)$ が魚 $u(x,t)$ を追いかける効果: $dv/dt$ に $ - \nabla ( v \nabla u)$ を足す

ことでそれぞれの効果をモデル化することができる．あとやるべきことはそれぞれの効果の「強さ」を例えば $\gamma_1$, $\gamma_2$ として調整してこれらの項にかけるぐらいかな．

すると、先同様にして $x \in [0, L] \times [0, L]$, $t > 0$ とし、以下のようなモデルを考えることになる．

$\left\{\begin{array}{rcl} \displaystyle \frac{\partial u}{\partial t} & = & \gamma_1 \nabla (u \nabla v) + \epsilon_1 \Delta u - C_1 uv + D_1 u, \cr\cr \displaystyle \frac{\partial v}{\partial t} & = & - \gamma_2 \nabla (v \nabla u) + \epsilon_2 \Delta v + C_2 uv - D_2 v \end{array}\right.$

ただし、前回と同じように、係数 $\gamma_1$, $\gamma_2$, $\epsilon_1$, $\epsilon_2$, $C_1$, $C_2$, $D_1$, $D_2$ は全て定数で正としておこう．

そして、空間領域に関する境界条件としては周期的境界条件を課しておくことにしよう．

method of line を適用する

ここでは、空間領域 $\Omega = [0, L] \times [0, L]$ に対し、その分割数を $N \times N$ とし、境界条件として周期的境界条件を離散化しよう．

よって、$\Delta x = \Delta y = L / N$ として $u(x,t)$ の離散化に相当する従属変数を $ u(t) = \{ u_{i,j}(t)\}_{i, j = 1}^N$ とし、境界条件は外側にはみ出す点を仮想的に考えたとして

\[ \left\{ \begin{array}{rcl} \displaystyle u_{0, j} & := & u_{N, j} {\rm \,\,\, for \,\,\,} \forall j ,\cr\cr \displaystyle u_{N+1, j} & := & u_{1, j} {\,\,\, 同上 } , \cr\cr \displaystyle u_{i, 0} & := & u_{i, N} {\rm \,\,\, for \,\,\,} \forall i, \cr\cr \displaystyle u_{i, N+1} & := & u_{i, 1} {\,\,\, 同上 } \end{array} \right. \]

と解釈することになる．

こうしておいて、まず、$\nabla (u \nabla v)$ の差分近似を考えてみよう．もちろんいろんなやり方があるが、たとえば、わかりやすいように、

\[ \nabla (u \nabla v) = \frac{\partial}{\partial x} \left( u \frac{\partial v}{\partial x} \right) + \frac{\partial}{\partial y} \left( u \frac{\partial v}{\partial y} \right) \]

の右辺第一項の離散化を

\[ \left\{ \left(\frac{u_{i+1,j} + u_{i,j}}{2} \right)\left(\frac{v_{i+1,j} - v_{i,j}}{\Delta x} \right) - \left(\frac{u_{i,j} + u_{i-1,j}}{2} \right)\left(\frac{v_{i,j} - v_{i-1,j}}{\Delta x} \right) \right\} \frac{1}{\Delta x} \]

とすると、これは等価変形して

\[ \frac{1}{2}\left\{ \left( \frac{u_{i+1,j} - u_{i,j} }{ \Delta x} \right)\left( \frac{ v_{i+1,j} - v_{i,j}}{\Delta x} \right) + \left( \frac{u_{i,j} - u_{i-1,j} }{ \Delta x} \right)\left( \frac{ v_{i,j} - v_{i-1,j}}{\Delta x} \right) \right\} + u\left( \frac{ v_{i+1,j} - 2v_{i,j} + v_{i-1,j} }{\Delta x^2} \right) \]

となる．あとは $y$ 成分も同様に計算して足そう．

すると、最初の中括弧の項の合計を($u,v$ について対称なことに留意して) $S_g(u,v)_{i,j}$ :

\[ \begin{array}{rcl} S_g(u,v)_{i,j} & := & \displaystyle \frac{1}{2}\left\{ \left( \frac{u_{i+1,j} - u_{i,j} }{ \Delta x} \right)\left( \frac{ v_{i+1,j} - v_{i,j}}{\Delta x} \right) + \left( \frac{u_{i,j} - u_{i-1,j} }{ \Delta x} \right)\left( \frac{ v_{i,j} - v_{i-1,j}}{\Delta x} \right)\right. \cr\cr % & & \displaystyle \left. + \left( \frac{u_{i,j+1} - u_{i,j} }{ \Delta y} \right)\left( \frac{ v_{i,j+1} - v_{i,j}}{\Delta y} \right) + \left( \frac{u_{i,j} - u_{i,j-1} }{ \Delta y} \right)\left( \frac{ v_{i,j} - v_{i,j-1}}{\Delta y} \right) \right\} \end{array} \]

と書くことにして、 $\nabla (u \nabla v)$ は

\[ \nabla (u \nabla v) \cong S_g(u,v)_{i,j} + u_{i,j} \left( \Delta v \right)_{i,j} \]

と離散近似できることがわかる． $\nabla (u \nabla v)$ も同様に、

\[ \nabla (v \nabla u) \cong S_g(u,v)_{i,j} + v_{i,j} \left( \Delta u \right)_{i,j} \]

と離散近似できる．

ちなみに、 $\left( \Delta u \right)_{i,j}$ は空間二階微分であるラプラシアンを $u$ に適用したものの離散近似で、前回紹介している．定義は

\[ (\Delta u)_{i,j} := \frac{ u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j+1} - 4 u_{i,j} }{\Delta x^2} \]

だった．

これで離散化すべき項の離散化が終わったので、これらを使った場合は、全体は $i = 1,2,\cdots, N$, $j = 1,2,\cdots, N$ に対して

\[ \left\{ \begin{array}{rcl} \displaystyle \frac{d u_{i,j}(t) }{dt} & = & \gamma_1 \left( S_g(u,v)_{i,j} + u_{i,j} \left( \Delta v \right)_{i,j} \right) + \epsilon_1 (\Delta u)_{i,j} - C_1 u_{i,j} v_{i,j} + D_1 u_{i,j}, \cr\cr \displaystyle \frac{d v_{i,j}(t) }{dt} & = & - \gamma_2 \left( S_g(u,v)_{i,j} + v_{i,j} \left( \Delta u \right)_{i,j} \right) + \epsilon_2 (\Delta v)_{i,j} + C_2 u_{i,j} v_{i,j} - D_2 v_{i,j}, \cr \end{array} \right. \]

という $N \times N$ 元連立常微分方程式で近似できることになる．

あとは少しずつプログラムを作っていくだけだ．まずは、問題のパラメータを設定してしまおう．

1
2
3
4
5
6
7
8


γ1, γ2 = 0.04, 0.01
ϵ1, ϵ2 = 0.001, 0.001
C1 = C2 = D1 = D2 = 1.0

L = 1.0
N = 40
Δx = Δy = L/N
Δt = 0.05 * Δx^2

なお、それっぽい挙動を数値計算結果でみるためのパラメータ調整は意外に難しい．試行錯誤が必要なので、自力で似たようなことをやるときはそれなりに時間と手間がかかるものと思っておこう．

次に、ラプラシアンを前回同様に構成しよう．それにはまず、係数行列を下記のように定義する．

1
2
3
4
5
6
7
8
9


using LinearAlgebra
using SparseArrays

d = [ -2.0 for i in 1:N ]
v = [ 1 for i in 1:N-1 ]
d2 = Diagonal(d) + diagm(1 => v, -1 => v)
d2[1,N] = d2[N, 1] = 1.0

d2 = sparse( (1/(Δx)^2) * d2 )

40×40 SparseMatrixCSC{Float64,Int64} with 120 stored entries: [1 , 1] = -3200.0 [2 , 1] = 1600.0 [40, 1] = 1600.0 [1 , 2] = 1600.0 [2 , 2] = -3200.0 [3 , 2] = 1600.0 [2 , 3] = 1600.0 [3 , 3] = -3200.0 [4 , 3] = 1600.0 [3 , 4] = 1600.0 [4 , 4] = -3200.0 [5 , 4] = 1600.0 ⋮ [36, 37] = 1600.0 [37, 37] = -3200.0 [38, 37] = 1600.0 [37, 38] = 1600.0 [38, 38] = -3200.0 [39, 38] = 1600.0 [38, 39] = 1600.0 [39, 39] = -3200.0 [40, 39] = 1600.0 [1 , 40] = 1600.0 [39, 40] = 1600.0 [40, 40] = -3200.0

こうしておいて、近似データそのものである行列 $u$ に対して、上で定義した行列 $d2$ を用いた $d2 * u + u * d2$ は $u$ に(離散)ラプラシアンを適用した結果になるというのは前回紹介したな．あとは使いやすいようにこれを関数にしておこう．

1

diff2(u) = u * d2 + d2 * u

次に、$S_g(u,v)_{i,j}$ を計算する関数を作ろう．少し長いが、次のようにすると良いだろう．丁寧に読めば、なにをやっているかはわかるはずだ．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45


# 周期的境界条件下で、u_{i+1,j} - u_{i,j} という感じの計算をする.
# 行列を渡して、行列が返ってくる．
# dim は行列の次元番号．1 で行方向(グラフの y方向)、
#                       2 で列方向(グラフの x方向)の計算になる．
function cyclic_diff(u, dim)
    res = diff(u, dims = dim)

    if dim == 1
        res = vcat(res, (u[1,:] - u[end,:])')
        # ベクトルになると強制的に縦ベクトルになってしまうので ' で戻している
    else
        res = hcat(res, u[:,1] - u[:,end])
    end
    return res
end    

# 周期的境界条件下でベクトルの添字を一つずらす
# 行列を渡して、行列が返ってくる．
# dim は行列の次元番号．1 で行方向(グラフの y方向)、
#                       2 で列方向(グラフの x方向)のずらしになる．
function cyclic_shift(u, dim)
    if dim == 1
        res = vcat(u[end,:]', u[1:end-1,:])
        # ベクトルになると強制的に縦ベクトルになってしまうので ' で戻している
    else
        res = hcat(u[:,end], u[:,1:end-1])
    end
    return res
end
            
# S_g(u,v)_{i,j} を計算する．
# 行列 2つを渡して、行列 1つが返ってくる．
# dim は行列の次元番号．1 で行方向(グラフの y方向)、
#                       2 で列方向(グラフの x方向)の計算になる．
function sqrgrad(u,v)
    res = zeros(size(u))
    coeff = 0.5 * [ 1/Δx^2, 1/Δy^2 ]
    
    for dim in 1:2
      dudv = cyclic_diff(u,dim) .* cyclic_diff(v,dim)
      res += coeff[dim] * ( dudv  + cyclic_shift(dudv, dim) )
    end

    return res
end

よって、微分方程式の右辺の近似計算式は次のようになるだろう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12


function f(U)
    u, v = U[:, 1:N], U[:, N+1:end]
    uv = u .* v
    d2u, d2v = diff2(u), diff2(v)
    sguv = sqrgrad(u,v)
    
    fu =   γ1 * ( sguv + u .* d2v ) + ϵ1 * d2u - C1 * uv + D1 * u
    fv = - γ2 * ( sguv + v .* d2u ) + ϵ2 * d2v + C2 * uv - D2 * v
    
    res = hcat( fu, fv )
    return res
end

ちなみに、上をよく見ると分かるように、今回は $u(x)$ の近似行列 $u_{i,j}$ と $v(x)$ の近似行列 $v_{i,j}$ を受け渡しする時は、左右につなげて、一つにして扱っている．上下でなく左右にしたのは、同時にグラフにするとき便利だから、だ．

さて、これで連立常微分方程式の右辺(の近似)が計算できたことになるので、いつものように Runge-Kutta 法を定義する．

1
2
3
4
5
6
7


function RK(u)
    f1 = f(u)
    f2 = f(u + Δt/2 * f1)
    f3 = f(u + Δt/2 * f2)
    f4 = f(u + Δt * f3)
    return u + Δt * (f1 + 2*f2 + 2*f3 + f4)/6
end

これで計算手順の準備は整った．

次に、前回同様の初期値を作ろう．以下のようになる．

1
2
3


u0 = [ 0.3+0.4*exp(-4.6*4*((i-N/2)^2 +(j-N/2)^2) / (N^2)) for i in 1:N, j in 1:N ]
v0 = [ 0.4+2.0*exp(-4.6*16*((i-N/4)^2 +(j-N/4)^2) / (N^2)) for i in 1:N, j in 1:N ]
U0 = hcat( u0, v0 )

40×80 Array{Float64,2}: 0.300099 0.300152 0.300227 0.300332 … 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300152 0.300232 0.300347 0.300507 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300227 0.300347 0.300519 0.300759 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300332 0.300507 0.300759 0.301109 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300474 0.300725 0.301084 0.301584 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300661 0.301012 0.301513 0.302211 … 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300902 0.30138 0.302064 0.303016 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.301202 0.301839 0.302751 0.304021 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.301566 0.302396 0.303584 0.305238 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.301994 0.303051 0.304563 0.306669 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.302481 0.303796 0.305677 0.308298 … 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.303016 0.304616 0.306903 0.310089 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.303584 0.305485 0.308203 0.311989 0.4 0.4 0.4 0.4 0.4 0.4 0.4 ⋮ ⋱ ⋮
0.302481 0.303796 0.305677 0.308298 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.301994 0.303051 0.304563 0.306669 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.301566 0.302396 0.303584 0.305238 … 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.301202 0.301839 0.302751 0.304021 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300902 0.30138 0.302064 0.303016 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300661 0.301012 0.301513 0.302211 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300474 0.300725 0.301084 0.301584 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300332 0.300507 0.300759 0.301109 … 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300227 0.300347 0.300519 0.300759 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300152 0.300232 0.300347 0.300507 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300099 0.300152 0.300227 0.300332 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.300063 0.300097 0.300145 0.300212 0.4 0.4 0.4 0.4 0.4 0.4 0.4

念の為、(既に知ってはいるが)初期値を目視しておこう．まずは u の初期値で、

1
2
3
4


using Plots

X = Y = Δx:Δx:L
wireframe(X, Y, u0)

png

次は v の初期値だ．

1

wireframe(X, Y, v0)

png

よし、前回と同じであることが確認できたな．

さて、これで準備は整ったので、早速、計算してみよう. ちなみに下記の数字だと，PC で2分以下，JuliaBox だと約 3分かかるだろう．授業時に待つのが大変という場合は，自分で適当に数字を調整しよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


using ProgressMeter

U = copy(U0)
sq_U = copy(U0)  

Nstep = 200000 # 全ステップ数
Nframe = 400   # とっておきたいデータ数
interval = Integer(floor(Nstep/Nframe))  # 何回おきに記録するか？を計算

@showprogress for i in 1:Nstep
    global U, sq_U
    U = RK(U)
    if i % interval == 0  # 一定数ステップおきに
        sq_U = cat(sq_U, U, dims = 3)  # データを 3次元目に追加していく
    end
end

Progress: 100%|█████████████████████████████████████████| Time: 0:03:05

試しに結果の数字をちょろっと見てみる．

1

sq_U[:, :, end]

40×80 Array{Float64,2}: 0.327783 0.320195 0.311539 0.303557 … 0.850092 0.857623 0.858945 0.320195 0.316362 0.312025 0.308787 0.848476 0.85247 0.850565 0.311539 0.312025 0.313075 0.315759 0.841477 0.841975 0.837476 0.303557 0.308787 0.315759 0.324705 0.830552 0.828244 0.822178 0.298531 0.308472 0.321034 0.335516 0.817949 0.813828 0.807231 0.298469 0.312462 0.329409 0.347818 … 0.805757 0.800723 0.794293 0.303884 0.320813 0.3405 0.360932 0.795031 0.7897 0.783743 0.313259 0.331962 0.352875 0.373738 0.785687 0.780505 0.77518 0.323989 0.343565 0.364622 0.384858 0.777309 0.772717 0.768235 0.333656 0.353537 0.374055 0.393015 0.769794 0.766264 0.762888 0.340454 0.360298 0.379863 0.397156 … 0.763153 0.761162 0.759155 0.343343 0.362884 0.381227 0.396647 0.75716 0.75715 0.756727 0.342341 0.361268 0.378172 0.391646 0.751455 0.753771 0.755008 ⋮ ⋱
0.22333 0.239967 0.260119 0.281609 0.760699 0.756771 0.757356 0.237918 0.246592 0.260778 0.278315 0.754264 0.751626 0.753864 0.25749 0.258173 0.265767 0.278334 … 0.752554 0.751561 0.755291 0.279645 0.273106 0.274056 0.280923 0.756035 0.75669 0.761538 0.300485 0.288098 0.283066 0.284179 0.763622 0.765827 0.771474 0.316933 0.300402 0.290608 0.286614 0.774017 0.777846 0.784213 0.327944 0.309112 0.29605 0.288046 0.786209 0.791947 0.799153 0.334238 0.314888 0.30003 0.289235 … 0.799448 0.807376 0.815482 0.337205 0.318822 0.303402 0.290934 0.813048 0.823107 0.831804 0.3379 0.321509 0.306477 0.293332 0.826177 0.837701 0.846173 0.336665 0.32286 0.309022 0.296237 0.837701 0.849413 0.856487 0.333352 0.322485 0.310708 0.299556 0.846173 0.856487 0.861029

あまり変なことにはなっていなさそうなので、プロットしてみよう．

まずは初期値を念のために…

1
2
3


range = (minimum(sq_U), maximum(sq_U))
X2 = Δx:Δx:2L
contour(X2, Y, sq_U[:, :, 1], clims=range, fill=true, aspect_ratio = 1)

png

ちなみに向かって左が魚の初期分布 $u(x,y,0)$ で、右がサメの初期分布 $v(x,y,0)$ になる．

次に、最終ステップ時の結果は次のようになるな．

1

contour(X2, Y, sq_U[:, :, end], clims=range, fill=true, aspect_ratio = 1)

png

ふむ、「すごく変な」ことにはなっていないようだな、とりあえず．

さて、前回同様、次にこれらを動画にして、わかりやすくみてみようではないか．

まず、たくさんの画像ファイルを入れておくディレクトリをコマンドで作らせよう．これは一回だけやれば良い．

1

run(`mkdir UV`)

次に、データを画像ファイルとして保存するコマンドを作って、それを使って保存しよう．まずは次のような感じでコマンドを作る．

1
2
3
4
5
6


using Printf # すぐ下の @sprintf を使いたいので．

function figureUV(squ, num)
  contour(X2, Y, squ[:, :, num], clims=range, fill=true, aspect_ratio = 1)
  savefig( "UV/" * @sprintf("%04d", num) * ".png" )
end

それから、実際に画像にして保存だ．

1
2
3


@showprogress for n in 1:size(sq_U)[3]
    figureUV(sq_U, n)
end

Progress: 100%|█████████████████████████████████████████| Time: 0:00:25

ffmpeg で動画に

次に、前回同様、ffmpeg で動画にする．以下のようにすればよいだろう．この過程は結構早いはず.

1

run(`ffmpeg -y -r 10 -i "UV/%04d.png" -pix_fmt yuv420p -s 1200x800 movie-uv.mp4`)

ffmpeg version N-80522-ge0faad8 Copyright © 2000-2016 the FFmpeg developers built with gcc 5.4.0 (GCC) configuration: –enable-gpl –enable-version3 –disable-w32threads –enable-dxva2 –enable-libmfx –enable-nvenc –enable-avisynth –enable-bzlib –enable-fontconfig –enable-frei0r –enable-gnutls –enable-iconv –enable-libass –enable-libbluray –enable-libbs2b –enable-libcaca –enable-libfreetype –enable-libgme –enable-libgsm –enable-libilbc –enable-libmodplug –enable-libmp3lame –enable-libopencore-amrnb –enable-libopencore-amrwb –enable-libopenjpeg –enable-libopus –enable-librtmp –enable-libschroedinger –enable-libsnappy –enable-libsoxr –enable-libspeex –enable-libtheora –enable-libtwolame –enable-libvidstab –enable-libvo-amrwbenc –enable-libvorbis –enable-libvpx –enable-libwavpack –enable-libwebp –enable-libx264 –enable-libx265 –enable-libxavs –enable-libxvid –enable-libzimg –enable-lzma –enable-decklink –enable-zlib libavutil 55. 24.100 / 55. 24.100 libavcodec 57. 46.100 / 57. 46.100 libavformat 57. 38.102 / 57. 38.102 libavdevice 57. 0.101 / 57. 0.101 libavfilter 6. 46.102 / 6. 46.102 libswscale 4. 1.100 / 4. 1.100 libswresample 2. 1.100 / 2. 1.100 libpostproc 54. 0.100 / 54. 0.100 Input #0, image2, from ‘UV/%04d.png’: Duration: 00:00:16.04, start: 0.000000, bitrate: N/A Stream #0:0: Video: png, rgb24(pc), 2700x1800, 25 fps, 25 tbr, 25 tbn, 25 tbc [libx264 @ 000000000262b980] using cpu capabilities: MMX2 SSE2Fast SSSE3 SSE4.2 [libx264 @ 000000000262b980] profile High, level 3.2 [libx264 @ 000000000262b980] 264 - core 148 r2694 3b70645 - H.264/MPEG-4 AVC codec - Copyleft 2003-2016 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=-2 threads=12 lookahead_threads=2 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=10 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=23.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00 [mp4 @ 0000000003d52320] Using AVStream.codec to pass codec parameters to muxers is deprecated, use AVStream.codecpar instead. Output #0, mp4, to ‘movie-uv.mp4’: Metadata: encoder : Lavf57.38.102 Stream #0:0: Video: h264 (libx264) ([33][0][0][0] / 0x0021), yuv420p, 1200x800, q=-1–1, 10 fps, 10240 tbn, 10 tbc Metadata: encoder : Lavc57.46.100 libx264 Side data: cpb: bitrate max/min/avg: 0/0/0 buffer size: 0 vbv_delay: -1 Stream mapping: Stream #0:0 -> #0:0 (png (native) -> h264 (libx264)) Press [q] to stop, [?] for help frame= 401 fps= 24 q=-1.0 Lsize= 1017kB time=00:00:39.80 bitrate= 209.3kbits/s speed=2.41x
video:1011kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 0.544998% [libx264 @ 000000000262b980] frame I:2 Avg QP:12.79 size: 12264 [libx264 @ 000000000262b980] frame P:101 Avg QP:16.28 size: 4343 [libx264 @ 000000000262b980] frame B:298 Avg QP:16.18 size: 1919 [libx264 @ 000000000262b980] consecutive B-frames: 0.7% 0.5% 0.0% 98.8% [libx264 @ 000000000262b980] mb I I16..4: 38.8% 56.1% 5.1% [libx264 @ 000000000262b980] mb P I16..4: 8.7% 11.1% 0.3% P16..4: 14.2% 3.9% 1.1% 0.0% 0.0% skip:60.9% [libx264 @ 000000000262b980] mb B I16..4: 1.4% 0.9% 0.0% B16..8: 16.9% 2.6% 0.1% direct: 1.3% skip:76.8% L0:48.1% L1:49.7% BI: 2.2% [libx264 @ 000000000262b980] 8x8 transform intra:51.2% inter:99.2% [libx264 @ 000000000262b980] coded y,uvDC,uvAC intra: 8.9% 39.1% 8.7% inter: 1.6% 4.3% 0.0% [libx264 @ 000000000262b980] i16 v,h,dc,p: 37% 40% 4% 19% [libx264 @ 000000000262b980] i8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 33% 26% 35% 1% 1% 1% 1% 1% 1% [libx264 @ 000000000262b980] i4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 41% 34% 19% 1% 2% 1% 1% 1% 1% [libx264 @ 000000000262b980] i8c dc,h,v,p: 46% 24% 20% 9% [libx264 @ 000000000262b980] Weighted P-Frames: Y:0.0% UV:0.0% [libx264 @ 000000000262b980] ref P L0: 64.9% 1.6% 22.5% 10.9% [libx264 @ 000000000262b980] ref B L0: 85.0% 12.2% 2.9% [libx264 @ 000000000262b980] ref B L1: 95.0% 5.0% [libx264 @ 000000000262b980] kb/s:206.47

これで movie-uv.mp4 という動画ファイルができたはずだ． julia/jyupiter のファイル操作画面で、その動画ファイルを view したり、ダウンロードして(jupyter はファイルを1つずつなら指定してから download できるぞ)動かしてみるなどして、確認しよう．

参考のために、そのファイルをここに直接おいておく．

得られる動画ファイル (mp4 形式)
得られる動画ファイル(色を強調した版) (mp4 形式)

ちなみに、contour の中の clims を手でいじることで色の濃さが変わることを利用したのが、上の二番目のファイルの「色強調版」だ．こちらのほうが見易いだろう．

そして、(とくに色強調版だと)魚がサメのいる方から逃げようとしていることや、「さざ波」のようなものが見えることがわかるだろう．

初期値を乱数で生成した場合

初期値作成部分を

1
2
3


u0 = rand(N,N)
v0 = rand(N,N)
U0 = hcat( u0, v0 )

として同様の「計算」「動画作成」をやってみると、以下のような感じだ．

計算過程等は上とまったく同じなので省略しよう．ちなみに、初期値を contour 表示すると

svg

というような感じになるはずだ．乱数で生成するから、人によって異なるがな．

そして、この初期値から動画にしたものもここに載せておこう．

得られる動画ファイル(初期値が乱数の場合) (mp4 形式)

途中の挙動の「もぞもぞした感じ」がなんとも生物っぽいような気がしないかな？

Report No.13 (チャレンジング: ちょっと難しい)

上のモデル偏微分方程式の中の常微分方程式部分、つまり、もともとの Volterra 方程式に存在した部分、この部分を改良して、上のように数値計算をし、動画を作ろう．ちなみにどのように改良するかであるが、wikipedia のロトカ・ヴォルテラの方程式を参照すると、「モデルの改良」という項目があるのでそこを参照するなどすればよいだろう．
上のストーリーを拡張して、3つの生物種の被食者・捕食者関係が循環するような状況を考えよう．つまり、生物として A, B, C の 3種類が居るとして、A は B を食べ、B は C を食べ、C は A を食べる、という状況を考える(ほんとにそんな状況があるかどうかは別として)．そしてその状況下で、上のストーリーのようなモデル方程式を作り、数値計算をして、挙動を動画化して様子をみてみよう．
なんでも良いので、なにかの問題に対して自分で偏微分方程式でモデルを作成して計算し、できれば動画を作成し、自分なりに納得できるか、評価してみよう．