格雷斯科特模型答案

【问题标题】：Gray Scott Model格雷斯科特模型
【发布时间】：2019-04-21 18:12:37
【问题描述】：

我试图展示扩散的灰色斯科特模型。即使我觉得我的代码非常接近正确，我仍然会收到一堆运行时警告错误。我的离散化有问题吗？

import numpy as np
import matplotlib.pyplot as plt

#parameters
N=128
F=.042
k=.062
Du=(2**-5)*(N**2/6.25)
Dv=(1**-5)*(N**2/6.25)
tend=1000                             
dt=tend/N
t=0

#start arrays
U=np.ones((N,N))
V=np.zeros((N,N))

#Initial Value Boxes (20x20 in middle)
low=int(((N/2)-10))
high=int(((N/2)+10))+1
U[low:high,low:high]=.5
V[low:high,low:high]=.25

#Random Noise
U+=.01*np.random.random((N,N))
V+=.01*np.random.random((N,N))

#Laplace
def Laplace(f):
    return np.roll(f,1)+np.roll(f,-1)+np.roll(f,1,axis=False)+np.roll(f,-1,axis=False)-4*f

#Solve
pstep=100
for t in range(tend):
    U+=((Du*Laplace(U))-(U*V*V)+(F*(1-U)))
    V+=((Dv*Laplace(V))+(U*V*V)-(V*(F+k)))
    if t%pstep ==0:
        print(t//pstep)
        plt.imshow(U, interpolation='bicubic',cmap=plt.cm.jet)

【问题讨论】：

一个问题是，您应该分别计算dU 和dV，然后分别计算U+=dU 和V+=dV，否则您使用更新后的U 来计算V。但这并不t 完全解决了数值不稳定性。我认为您可能计算 Du 和 Dv 错误？价值似乎很高

标签： python numpy pde

【解决方案1】：

好的，我通过改变计算中的一些东西来让它工作，但主要是通过大量降低扩散系数和减少时间步长来改变数值稳定性。这样做的最终结果是，整个模拟在每一步之间变化较小，因此变化的值要小得多。

您遇到的错误是由于计算 dU 和 dV 时的浮点数溢出，通过减慢整个过程（更多时间步）您不需要在 dU 和dV

import numpy as np
import matplotlib.pyplot as plt

# parameters
N = 128
F = .042
k = .062
# Du=(2**-5)*(N**2/6.25) # These were way too high for the 
# numeric stability given the timestep
Du = 0.1
# Dv=(1**-5)*(N**2/6.25)
Dv = 0.5
tend = 1000
dt = tend / N
t = 0
dt = 0.1  # Timestep - 
# the smaller you go here, the faster you can let the diffusion coefficients be

# start arrays
U = np.ones((N, N))
V = np.zeros((N, N))

# Initial Value Boxes (20x20 in middle)
low = int(((N / 2) - 10))
high = int(((N / 2) + 10)) + 1
U[low:high, low:high] = .5
V[low:high, low:high] = .25

# Random Noise
U += .01 * np.random.random((N, N))
V += .01 * np.random.random((N, N))


# Laplace
def Laplace(f):
    return np.roll(f, 1) + np.roll(f, -1) + np.roll(f, 1, axis=False) + np.roll(f,-1,                                                                              axis=False) - 4 * f


# Solve
pstep = 100
for t in range(tend):
    dU = ((Du * Laplace(U)) - (U * V * V) + (F * (1 - U))) * dt
    dV = ((Dv * Laplace(V)) + (U * V * V) - (V * (F + k))) * dt
    U += dU
    V += dV
    if t % pstep == 0:
        print(t // pstep)
        plt.imshow(U, interpolation='bicubic', cmap=plt.cm.jet, vmin=0, vmax=1)

当然，我所做的更改会稍微改变物理特性，您需要更改您的 t 和 pstep 以使这些值有意义。还要检查你是如何计算Du 和Dv 的。如果这些值实际上应该是 ~8000，那么你需要一个小得多的时间步长。

供其他人参考，格雷斯科特模型解释为here

【讨论】：