这是一个速度问题。我试图求解一个扩散方程，它有三种行为状态，其中：
- Lambda==0平衡
- λ>0最大扩散
- Lambda<0分钟扩散
瓶颈是扩散算子函数中的其他语句
平衡态有一个简单的T算子和扩散算子。对于其他两个州来说，情况变得相当复杂。到目前为止，我还没有耐心看完代码运行时间。据我所知，方程是正确的，平衡态的输出似乎是正确的。也许有人有一些技巧可以提高非平衡态的速度？
（欧拉解（FTCS）而不是龙格-库塔会更快，我想。还没有尝试过。）
你可以导入任何黑白图像来尝试代码。
import numpy as np
import sympy as sp
import scipy.ndimage.filters as flt
from PIL import Image
# import image
im = Image.open("/home/will/Downloads/zebra.png")
arr = np.array(im)
arr=arr/253.
def T(lamda,x):
"""
T Operator
lambda is a "steering" constant between 3 behavior states
-----------------------------
0 -> linearity
+inf -> max
-inf -> min
-----------------------------
"""
if lamda == 0: # linearity
return x
elif lamda > 0: # Half-wave rectification
return np.max(x,0)
elif lamda < 0: # Inverse half-wave rectification
return np.min(x,0)
def Diffusion_operator(lamda,f,t):
"""
2D Spatially Discrete Non-Linear Diffusion
------------------------------------------
Special case where lambda == 0, operator becomes Laplacian
Parameters
----------
D : int diffusion coefficient
h : int step size
t0 : int stimulus injection point
stimulus : array-like luminance distribution
Returns
----------
f : array-like output of diffusion equation
-----------------------------
0 -> linearity (T[0])
+inf -> positive K(lamda)
-inf -> negative K(lamda)
-----------------------------
"""
if lamda == 0: # linearity
return flt.laplace(f)
else: # non-linearity
f_new = np.zeros(f.shape)
for x in np.arange(0,f.shape[0]-1):
for y in np.arange(0,f.shape[1]-1):
f_new[x,y]=T(lamda,f[x+1,y]-f[x,y]) + T(lamda,f[x-1,y]-f[x,y]) + T(lamda,f[x,y+1]-f[x,y])
+ T(lamda,f[x,y-1]-f[x,y])
return f_new
def Dirac_delta_test(tester):
# Limit injection to unitary multiplication, not infinite
if np.sum(sp.DiracDelta(tester)) == 0:
return 0
else:
return 1
def Runge_Kutta(stimulus,lamda,t0,h,N,D,t_N):
"""
4th order Runge-Kutta solution to:
linear and spatially discrete diffusion equation (ignoring leakage currents)
Adiabatic boundary conditions prevent flux exchange over domain boundaries
Parameters
---------------
stimulus : array_like input stimuli [t,x,y]
lamda : int 0 +/- inf
t0 : int point of stimulus "injection"
h : int step size
N : int array size (stimulus.shape[1])
D : int diffusion coefficient [constant]
Returns
----------------
f : array_like computed diffused array
"""
f = np.zeros((t_N+1,N,N)) #[time, equal shape space dimension]
t = np.zeros(t_N+1)
if lamda ==0:
""" Linearity """
for n in np.arange(0,t_N):
k1 = D*flt.laplace(f[t[n],:,:]) + stimulus*Dirac_delta_test(t[n]-t0)
k1 = k1.astype(np.float64)
k2 = D*flt.laplace(f[t[n],:,:]+(0.5*h*k1)) + stimulus*Dirac_delta_test((t[n]+(0.5*h))- t0)
k2 = k2.astype(np.float64)
k3 = D*flt.laplace(f[t[n],:,:]+(0.5*h*k2)) + stimulus*Dirac_delta_test((t[n]+(0.5*h))-t0)
k3 = k3.astype(np.float64)
k4 = D*flt.laplace(f[t[n],:,:]+(h*k3)) + stimulus*Dirac_delta_test((t[n]+h)-t0)
k4 = k4.astype(np.float64)
f[n+1,:,:] = f[n,:,:] + (h/6.) * (k1 + 2.*k2 + 2.*k3 + k4)
t[n+1] = t[n] + h
return f
else:
""" Non-Linearity THIS IS SLOW """
for n in np.arange(0,t_N):
k1 = D*Diffusion_operator(lamda,f[t[n],:,:],t[n]) + stimulus*Dirac_delta_test(t[n]-t0)
k1 = k1.astype(np.float64)
k2 = D*Diffusion_operator(lamda,(f[t[n],:,:]+(0.5*h*k1)),t[n]) + stimulus*Dirac_delta_test((t[n]+(0.5*h))- t0)
k2 = k2.astype(np.float64)
k3 = D*Diffusion_operator(lamda,(f[t[n],:,:]+(0.5*h*k2)),t[n]) + stimulus*Dirac_delta_test((t[n]+(0.5*h))-t0)
k3 = k3.astype(np.float64)
k4 = D*Diffusion_operator(lamda,(f[t[n],:,:]+(h*k3)),t[n]) + stimulus*Dirac_delta_test((t[n]+h)-t0)
k4 = k4.astype(np.float64)
f[n+1,:,:] = f[n,:,:] + (h/6.) * (k1 + 2.*k2 + 2.*k3 + k4)
t[n+1] = t[n] + h
return f
# Code to run
N=arr.shape[1] # Image size
stimulus=arr[0:N,0:N,1]
D = 0.3 # Diffusion coefficient [0>D>1]
h = 1 # Runge-Kutta step size [h > 0]
t0 = 0 # Injection time
t_N = 100 # End time
f_out_equil = Runge_Kutta(stimulus,0,t0,h,N,D,t_N)
f_out_min = Runge_Kutta(stimulus,-1,t0,h,N,D,t_N)
f_out_max = Runge_Kutta(stimulus,1,t0,h,N,D,t_N)
简言之，f_out_equil的计算速度相对较快，而最小和最大情况则昂贵且缓慢。
以下是我一直在使用的图像链接：http://4.bp.blogspot.com/_KbtOtXslVZE/SweZiZWllzI/AAAAAAAAAIg/i9wc-yfdW78/s200/Zebra_Black_and_White_by_Jenvanw.jpg
感谢关于改进我的编码的提示，非常感谢，
以下是输出的快速绘图脚本
import matplotlib.pyplot as plt
fig1, (ax1,ax2,ax3,ax4,ax5) = plt.subplots(ncols=5, figsize=(15,5))
ax1.imshow(f_out_equil[1,:,:],cmap='gray')
ax2.imshow(f_out_equil[t_N/10,:,:],cmap='gray')
ax3.imshow(f_out_equil[t_N/2,:,:],cmap='gray')
ax4.imshow(f_out_equil[t_N/1.5,:,:],cmap='gray')
ax5.imshow(f_out_equil[t_N,:,:],cmap='gray')
对于python中的循环往往很慢；你可以通过尽可能多的矢量化来获得巨大的加速。（这将对你今后的任何数值问题都有很大帮助）。新的T算子一次处理整个阵列，并且在Diffusion_operator中对np.roll的调用使图像阵列正确排列，以进行有限差分计算。
整个过程在我的电脑上运行了大约10秒。
def T(lamda,x):
"""
T Operator
lambda is a "steering" constant between 3 behavior states
-----------------------------
0 -> linearity
+inf -> max
-inf -> min
-----------------------------
"""
if lamda == 0: # linearity
return x
elif lamda > 0: # Half-wave rectification
maxval = np.zeros_like(x)
return np.array([x, maxval]).max(axis=0)
elif lamda < 0: # Inverse half-wave rectification
minval = np.zeros_like(x)
return np.array([x, minval]).min(axis=0)
def Diffusion_operator(lamda,f,t):
"""
2D Spatially Discrete Non-Linear Diffusion
------------------------------------------
Special case where lambda == 0, operator becomes Laplacian
Parameters
----------
D : int diffusion coefficient
h : int step size
t0 : int stimulus injection point
stimulus : array-like luminance distribution
Returns
----------
f : array-like output of diffusion equation
-----------------------------
0 -> linearity (T[0])
+inf -> positive K(lamda)
-inf -> negative K(lamda)
-----------------------------
"""
if lamda == 0: # linearity
return flt.laplace(f)
else: # non-linearity
f_new = T(lamda,np.roll(f,1, axis=0) - f)
+ T(lamda,np.roll(f,-1, axis=0) - f)
+ T(lamda,np.roll(f, 1, axis=1) - f)
+ T(lamda,np.roll(f,-1, axis=1) - f)
return f_new