Julia 挑战 - FitzHugh–Nagumo 模型 PDE Runge-Kutta 求解器
Julia challenge - FitzHugh–Nagumo model PDE Runge-Kutta solver
我是Julia编程语言的新手,所以我不太了解如何优化代码。我听说 Julia 应该比 Python 快,但是我写了一个简单的 Julia code for solving the FitzHugh–Nagumo model ,它似乎并不比 Python 快。
FitzHugh–Nagumo 模型方程为:
function FHN_equation(u,v,a0,a1,d,eps,dx)
u_t = u - u.^3 - v + laplacian(u,dx)
v_t = eps.*(u - a1 * v - a0) + d*laplacian(v,dx)
return u_t, v_t
end
其中u
和v
是变量,是2D字段(即二维数组),a0,a1,d,eps
是模型的参数。参数和变量都是 Float 类型。 dx
为控制网格点间距的参数,用于拉普拉斯函数的使用,拉普拉斯函数是周期性边界条件的有限差分的实现。
如果你们中的一位 Julia 编码专家能给我一些提示,告诉我如何在 Julia 中做得更好,我会很高兴听到。
Runge-Kutte 阶跃函数为:
function uv_rk4_step(Vs,Ps, dt)
u = Vs.u
v = Vs.v
a0=Ps.a0
a1=Ps.a1
d=Ps.d
eps=Ps.eps
dx=Ps.dx
du_k1, dv_k1 = FHN_equation(u,v,a0,a1,d,eps,dx)
u_k1 = dt*du_k1י
v_k1 = dt*dv_k1
du_k2, dv_k2 = FHN_equation((u+(1/2)*u_k1),(v+(1/2)*v_k1),a0,a1,d,eps,dx)
u_k2 = dt*du_k2
v_k2 = dt*dv_k2
du_k3, dv_k3 = FHN_equation((u+(1/2)*u_k2),(v+(1/2)*v_k2),a0,a1,d,eps,dx)
u_k3 = dt*du_k3
v_k3 = dt*dv_k3
du_k4, dv_k4 = FHN_equation((u+u_k3),(v+v_k3),a0,a1,d,eps,dx)
u_k4 = dt*du_k4
v_k4 = dt*dv_k4
u_next = u+(1/6)*u_k1+(1/3)*u_k2+(1/3)*u_k3+(1/6)*u_k4
v_next = v+(1/6)*v_k1+(1/3)*v_k2+(1/3)*v_k3+(1/6)*v_k4
return u_next, v_next
end
我使用 PyPlot 包中的 imshow() 来绘制 u 字段。
这不是一个完整的答案,而是对 laplacian
函数的优化尝试的尝试。 10x10 矩阵上的原始 laplacian
给了我 @time:
0.000038 seconds (51 allocations: 12.531 KB)
虽然这个版本:
function laplacian2(a,dx)
# Computes Laplacian of a matrix
# Usage: al=laplacian(a,dx)
# where dx is the grid interval
ns=size(a,1)
ns != size(a,2) && error("Input matrix must be square")
aa=zeros(ns+2,ns+2)
for i=1:ns
aa[i+1,1]=a[i,end]
aa[i+1,end]=a[i,1]
aa[1,i+1]=a[end,i]
aa[end,i+1]=a[1,i]
end
for i=1:ns,j=1:ns
aa[i+1,j+1]=a[i,j]
end
lap = Array{eltype(a),2}(ns,ns)
scale = inv(dx*dx)
for i=1:ns,j=1:ns
lap[i,j]=(aa[i,j+1]+aa[i+2,j+1]+aa[i+1,j]+aa[i+1,j+2]-4*aa[i+1,j+1])*scale
end
return lap
end
给@时间:
0.000010 seconds (6 allocations: 2.250 KB)
请注意分配减少。额外分配通常表示优化的潜力。
我是Julia编程语言的新手,所以我不太了解如何优化代码。我听说 Julia 应该比 Python 快,但是我写了一个简单的 Julia code for solving the FitzHugh–Nagumo model ,它似乎并不比 Python 快。
FitzHugh–Nagumo 模型方程为:
function FHN_equation(u,v,a0,a1,d,eps,dx)
u_t = u - u.^3 - v + laplacian(u,dx)
v_t = eps.*(u - a1 * v - a0) + d*laplacian(v,dx)
return u_t, v_t
end
其中u
和v
是变量,是2D字段(即二维数组),a0,a1,d,eps
是模型的参数。参数和变量都是 Float 类型。 dx
为控制网格点间距的参数,用于拉普拉斯函数的使用,拉普拉斯函数是周期性边界条件的有限差分的实现。
如果你们中的一位 Julia 编码专家能给我一些提示,告诉我如何在 Julia 中做得更好,我会很高兴听到。
Runge-Kutte 阶跃函数为:
function uv_rk4_step(Vs,Ps, dt)
u = Vs.u
v = Vs.v
a0=Ps.a0
a1=Ps.a1
d=Ps.d
eps=Ps.eps
dx=Ps.dx
du_k1, dv_k1 = FHN_equation(u,v,a0,a1,d,eps,dx)
u_k1 = dt*du_k1י
v_k1 = dt*dv_k1
du_k2, dv_k2 = FHN_equation((u+(1/2)*u_k1),(v+(1/2)*v_k1),a0,a1,d,eps,dx)
u_k2 = dt*du_k2
v_k2 = dt*dv_k2
du_k3, dv_k3 = FHN_equation((u+(1/2)*u_k2),(v+(1/2)*v_k2),a0,a1,d,eps,dx)
u_k3 = dt*du_k3
v_k3 = dt*dv_k3
du_k4, dv_k4 = FHN_equation((u+u_k3),(v+v_k3),a0,a1,d,eps,dx)
u_k4 = dt*du_k4
v_k4 = dt*dv_k4
u_next = u+(1/6)*u_k1+(1/3)*u_k2+(1/3)*u_k3+(1/6)*u_k4
v_next = v+(1/6)*v_k1+(1/3)*v_k2+(1/3)*v_k3+(1/6)*v_k4
return u_next, v_next
end
我使用 PyPlot 包中的 imshow() 来绘制 u 字段。
这不是一个完整的答案,而是对 laplacian
函数的优化尝试的尝试。 10x10 矩阵上的原始 laplacian
给了我 @time:
0.000038 seconds (51 allocations: 12.531 KB)
虽然这个版本:
function laplacian2(a,dx)
# Computes Laplacian of a matrix
# Usage: al=laplacian(a,dx)
# where dx is the grid interval
ns=size(a,1)
ns != size(a,2) && error("Input matrix must be square")
aa=zeros(ns+2,ns+2)
for i=1:ns
aa[i+1,1]=a[i,end]
aa[i+1,end]=a[i,1]
aa[1,i+1]=a[end,i]
aa[end,i+1]=a[1,i]
end
for i=1:ns,j=1:ns
aa[i+1,j+1]=a[i,j]
end
lap = Array{eltype(a),2}(ns,ns)
scale = inv(dx*dx)
for i=1:ns,j=1:ns
lap[i,j]=(aa[i,j+1]+aa[i+2,j+1]+aa[i+1,j]+aa[i+1,j+2]-4*aa[i+1,j+1])*scale
end
return lap
end
给@时间:
0.000010 seconds (6 allocations: 2.250 KB)
请注意分配减少。额外分配通常表示优化的潜力。