參考網路來源
https://ipython-books.github.io/124-simulating-a-partial-differential-equation-reaction-diffusion-systems-and-turing-patterns/
我想要把裡面的方程式，換成以下論文裡面的方程式
http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0174946
所以我目前的進度是這樣:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
mu = 0.16
nu= 0.04
D_H = 0.3
D_A=0.02
D_S=0.06
RHO_A = 0.03
RHO_H = 0.0001
epsilon=1.0
gamma=0.02
C=0.002
D=0.008
E=0.1
F=10
C_O=0.02
size = 50 # size of the 2D grid
dx = 2. / size # space step
T = 9.0 # total time
dt = .001 # time step
n = int(T / dt) # number of iterations
A = np.random.rand(size, size)
H = np.random.rand(size, size)
S = np.random.rand(size, size)
Y = np.random.rand(size, size)
def laplacian( Z, R, i , j , A ):
Z= 0.0 , dx=DELtA_X , dy=DELTA_X * 0.8 ;
if ( i > L ):
Z+= ( A[i-1][j] - A[i][j] )/dx/dx;
if ( i < R ):
Z+= ( A[i+1][j] - A[i][j] )/dx/dx;
if ( j > 0 ):
Z+= ( A[i][j-1] - A[i][j] )/dy/dy;
if ( j < Z ):
Z+= ( A[i][j+1] - A[i][j] )/dy/dy;
return Z;
#def laplacian(Z):
# Ztop = Z[0:-2, 1:-1]
# Zleft = Z[1:-1, 0:-2]
# Zbottom = Z[2:, 1:-1]
# Zright = Z[1:-1, 2:]
# Zcenter = Z[1:-1, 1:-1]
# return (Ztop + Zleft + Zbottom + Zright -
# 4 * Zcenter) / dx**2
def show_patterns(S, ax=None):
ax.imshow(S, cmap=plt.cm.copper,
interpolation='bilinear',
extent=[-1, 1, -1, 1])
ax.set_axis_off()
fig, axes = plt.subplots(3, 3, figsize=(8, 8))
step_plot = n // 9
# We simulate the PDE with the finite difference
# method.
for i in range(n):
# We compute the Laplacian of S and Y.
deltaA = laplacian(A)
deltaH = laplacian(H)
deltaS = laplacian(S)
# We take the values of S and Y inside the grid.
Ac = A[1:-1, 1:-1]
Hc = H[1:-1, 1:-1]
Sc = S[1:-1, 1:-1]
Yc = Y[1:-1, 1:-1]
# We update the variables.
A[1:-1, 1:-1],H[1:-1, 1:-1],S[1:-1, 1:-1], Y[1:-1, 1:-1]=\
Ac+dt * (C*A*A* S /H - mu * A + RHO_A * Y+D_A*deltaA),\
Hc+dt * (C * A*S - nu*H+ +RHO_H * Y+D_H*deltaH),\
Sc+dt * ( C_O - gamma * S - epsilon * Y * S+D_S*deltaS),\
Yc+dt * (D * A - E * Y + Y * Y / ( 1 + F * Y * Y))
# Neumann conditions: derivatives at the edges
# are null.
for Z in (Y, S):
if ( i > L ):
Z+= ( A[i-1][j] - A[i][j] )/dx/dx;
if ( i < R ):
Z+= ( A[i+1][j] - A[i][j] )/dx/dx;
if ( j > 0 ):
Z+= ( A[i][j-1] - A[i][j] )/dy/dy;
if ( j < Z ):
Z+= ( A[i][j+1] - A[i][j] )/dy/dy;
return Z;
# We plot the state of the system at
# 9 different times.
if i % step_plot == 0 and i < 9 * step_plot:
ax = axes.flat[i // step_plot]
show_patterns(S, ax=ax)
ax.set_title(f'$t={i * dt:.2f}$')
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
show_patterns(S, ax=ax)
可是laplacian出現問題，矩陣怪怪的。
論文裡面是用cuda的軟體寫的(我實在不知道怎麼改) :
__device__ float laplace_2( mat_t A, int i, int j , int L, int R, int Z){
float t = 0.0f, dx = DELTA_X, dy = (DELTA_X*0.8);
if ( i > L ) t+= ( A[i-1][j] - A[i][j] )/dx/dx;
if ( i < R ) t+= ( A[i+1][j] - A[i][j] )/dx/dx;
if ( j > 0 ) t+= ( A[i][j-1] - A[i][j] )/dy/dy;
if ( j < Z ) t+= ( A[i][j+1] - A[i][j] )/dy/dy;
return t;
}
最後我需要看S-Y curve，所以我又加了
#plot that shit
plt.figure()
plt.plot(V,U, label ='one' )
plt.xlabel('U')
plt.ylabel('V')
plt.legend()
plt.show()
可是又跟論文裡面 Fig1的那個圖案差很多。
我分離不出來單獨一條 V,U，所以是一大坨線，我希望畫單獨一條。