Julia Herausforderung - FitzHugh-Nagumo Modell PDE Runge-Kutta Löser

Question

Ich bin ein Neuling in der Programmiersprache Julia, daher weiß ich nicht viel darüber, wie man einen Code optimieren kann. Ich habe gehört, dass Julia im Vergleich zu Python schneller sein sollte, aber ich habe eine einfache geschrieben Julia code for solving the FitzHugh–Nagumo model, und es scheint nicht schneller als Python zu sein.

Die FitzHugh-Nagumo Modellgleichungen sind:

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 

wo uv und die Variablen, die 2D-Felder sind (das heißt, 2-dimensionale Arrays) und a0,a1,d,eps sind die Parameter des Modells. Beide Parameter und die Variablen sind vom Typ Float. dx ist der Parameter, der die Trennung zwischen Gitterpunkten für die Verwendung der Laplace-Funktion steuert, die eine Implementierung von finiten Differenzen mit periodischen Randbedingungen darstellt.

Wenn einer von Ihnen Experte Julia Coders mir einen Hinweis geben kann, wie man Dinge in Julia besser macht, werde ich mich freuen zu hören.

Die Runge-Kutte Sprungfunktion ist:

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 

Und ich habe verwendet imshow() aus PyPlot Paket, um das u Feld zu zeichnen.

Answer 1

Dies ist keine vollständige Antwort, sondern ein Vorgeschmack auf einen Optimierungsversuch auf die laplacian Funktion. Die ursprüngliche laplacian auf einer 10x10 Matrix gab mir die @time:

0.000038 seconds (51 allocations: 12.531 KB) 

Während dieser Version:

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 

Gibt @time:

0.000010 seconds (6 allocations: 2.250 KB) 

Beachten Sie die Reduzierung der Zuweisungen. Zusätzliche Zuordnungen zeigen normalerweise das Optimierungspotenzial an.