Q&A
def laplacian(ix, iy, s): の各入力変数の簡単な具体例と想定される出力結果を教えてください。
コードを読みましたがよく分かりません。
前提・実現したいこと
拡散方程式を用いてサバンナの植生を差分法でシミュレーションするコードを書いています．
発生している問題・エラーメッセージ
途中，配列を乗算する箇所で計算結果が欠損値(nan)になってしまいます．
計算はlaplacian()，calc_test()の部分で行い，問題もそこで発生しています．
以下に一例として，ラプラシアンの結果を示します．(ln)
[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.01 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0.01 -0.04 0.01 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0.01 0. 0. ] [ 0. 0. 0. 0.01 0. 0. 0. 0. 0. 0. ] [ 0. 0.01 0.01 -0.04 0.01 0. 0. 0. 0.01 0. ] [ 0. -0.04 0.01 0.01 0. 0. 0. 0.01 -0.04 0. ] [ 0. 0.01 0. 0. 0. 0.01 0. 0. 0.01 0. ] [ 0. 0. 0. 0. 0.01 -0.04 0.01 0.01 0. 0. ] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]] gs_00.png [[ 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00] [ 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 4.00000e-06 9.96480e-03 4.00000e-06 0.00000e+00] [ 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 2.00000e-06 9.96480e-03 -3.98832e-02 9.96480e-03 0.00000e+00] [ 0.00000e+00 0.00000e+00 0.00000e+00 2.00000e-06 0.00000e+00 0.00000e+00 4.00000e-06 9.96480e-03 4.00000e-06 0.00000e+00] [ 0.00000e+00 2.00000e-06 4.00000e-06 9.96480e-03 4.00000e-06 0.00000e+00 0.00000e+00 2.00000e-06 2.00000e-06 0.00000e+00] [ 0.00000e+00 9.96680e-03 9.96880e-03 -3.98832e-02 9.96480e-03 2.00000e-06 0.00000e+00 4.00000e-06 9.96480e-03 0.00000e+00] [ 0.00000e+00 -3.98852e-02 9.96880e-03 9.96680e-03 4.00000e-06 2.00000e-06 2.00000e-06 9.96480e-03 -3.98852e-02 0.00000e+00] [ 0.00000e+00 9.96480e-03 4.00000e-06 2.00000e-06 ... truncated gs_00.png [[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] gs_00.png [[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] gs_00.png [[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
該当のソースコード
Python
1import numpy as np 2import time 3from numba import jit 4from functools import wraps 5import matplotlib.pyplot as plt 6import matplotlib.animation as animation 7from matplotlib.ticker import MultipleLocator 8import random 9import sys 10 11@jit 12def laplacian(ix, iy, s): 13 ts = 0.0 14 ts += s[ix-1, iy] 15 ts += s[ix, iy-1] 16 ts += s[ix+1, iy] 17 ts += s[ix, iy+1] 18 ts -= s[ix, iy]*4.0 19 20 return ts 21 22@jit 23def calc_test(w, s, n, w2, s2, n2): 24 L = w.shape[0] 25 D_w = 0.01 26 D_s = 0.005 27 D_n = 0.001 28 dt = 0.2 29 f = 0.625 30 p = 0.05 31 q = 0.09004 32 m = 0.0096 33 A_r = 0.0 34 lw = np.zeros((L, L)) 35 ls = np.zeros((L, L)) 36 ln = np.zeros((L, L)) 37 w_x = np.zeros((L, L)) 38 for ix in range(1, L - 1): 39 for iy in range(1, L - 1): 40 lw[ix, iy] = D_w*laplacian(ix, iy, w) 41 ls[ix, iy] = D_s*laplacian(ix ,iy, s) 42 ln[ix, iy] = D_n*laplacian(ix, iy, n) 43 w_x[ix, iy] = A_r*differential_x(ix, iy, w) 44 45 cw = 1.0 - w - n*w*(1.0 - s/n/f) + w_x 46 cs = n*w*(1.0 - s/n/f) -p*n*s 47 cn = q*p*n*s - m*n 48 print(ln) 49 w2[:] = w + (lw + cw)*dt 50 s2[:] = s + (ls + cs)*dt 51 n2[:] = n + (ln + cn) * dt 52 53def runnning(): 54 steps = 5 55 stepbstep = steps/20 56 w2 = np.zeros((L, L)) 57 s2 = np.zeros((L, L)) 58 n2 = np.zeros((L, L)) 59 60 for i in range(steps): 61 if i%2 == 0: 62 calc_test(w, s, n, w2, s2, n2) 63 else: 64 calc_test(w2, s2, n2, w, s, n) 65 if i%stepbstep == 0: 66 filename = "gs_{:02d}.png".format(i//100) 67 print(filename) 68 plt.imshow(n, cmap='inferno') 69 plt.savefig(filename) 70 71L = 10 72w = np.zeros((L, L)) 73s = np.zeros((L, L)) 74n = np.zeros((L, L)) 75randomd_distribution(w, s, n) 76#center_distribution(w,5) 77#center_distribution(s,4) 78#center_distribution(n) 79runnning()
試したこと
ここに問題に対して試したことを記載してください。
補足情報（FW/ツールのバージョンなど）
Python 3.6.4
randomd_distributionが未定義です。実行可能なコードを提示ください。
あなたの回答
tips
プレビュー
