Julia 挑战 - FitzHugh–Nagumo 模型 PDE Runge-Kutta 求解器答案

【问题标题】：Julia challenge - FitzHugh–Nagumo model PDE Runge-Kutta solverJulia 挑战 - FitzHugh–Nagumo 模型 PDE Runge-Kutta 求解器
【发布时间】：2016-11-02 03:49:54
【问题描述】：

我是 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是变量，是二维字段（即二维数组），a0,a1,d,eps是模型的参数。参数和变量都是浮点类型。 dx是控制网格点间距的参数，用于拉普拉斯函数，它是具有周期性边界条件的有限差分的实现。

如果你们中的一位 Julia 编码专家可以给我一个提示，告诉我如何在 Julia 中做得更好，我会很高兴听到。

龙格-库特阶跃函数为：

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 字段。

【问题讨论】：

乍一看你应该删除向量化操作并将它们写成显式循环。
虽然对于速度目的来说不是必需的（假设代码格式正确），但对于这样的问题，在函数参数中包含类型信息通常很有用。这使得不熟悉 FitzHugh-Nagumo 的人更更容易就您的代码提供建议。
更好的是提供一个独立的子程序来模拟输入并计算您的函数时间，我们可以将其复制并粘贴到 REPL 中，然后尝试改进......这些是通常的技巧从而从 StackOverflow 中获得最佳答案。
支持@ColinTBowers cmets，几乎总是值得添加有关向量、矩阵的维度（不是元素类型）的类型信息。因为对于具有其他接口的类型，这些功能不会真正“工作”。这也为功能的人类读者提供了更好的界面。类型可以尽可能抽象，以保持功能尽可能通用。
最后的评论可能会节省您的时间。我愿意大打赌你来自 Matlab 背景（和我一样）。您链接的 github 似乎每个文件都有一个功能（或每个文件一个类型）。经典的 Matlab。对于 Julia（或大多数其他语言），您无需担心这一点。您可以将多个函数和类型放在一个文件中并使其成为一个模块。一旦你习惯了它方式会更方便。

标签： julia pde runge-kutta

【解决方案1】：

这不是一个完整的答案，而是对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

给@time：

0.000010 seconds (6 allocations: 2.250 KB)

注意分配的减少。额外的分配通常表明有优化的潜力。

【讨论】：