Julia challange - modèle FitzHugh-Nagumo solveur PDE Runge-Kutta

Question

Je suis débutant dans le langage de programmation Julia, donc je ne sais pas comment optimiser un code. J'ai entendu dire que Julia devrait être plus rapide par rapport à Python, mais j'ai écrit un simple Julia code for solving the FitzHugh–Nagumo model, et il ne semble pas être plus rapide que Python.

les équations du modèle FitzHugh-Nagumo sont:

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 

où u et v sont des variables, qui sont des champs 2D (à savoir, 2 tableaux à une dimension), et a0,a1,d,eps sont les paramètres du modèle. Les deux paramètres et les variables sont de type Float. dx est le paramètre qui contrôle la séparation entre le point de grille, pour l'utilisation de la fonction laplacienne, qui est une implémentation de différences finies avec des conditions aux limites périodiques.

Si l'un de vos codeurs experts Julia peut me donner un indice sur la façon de mieux faire les choses à Julia, je serai heureux de l'entendre.

La fonction étape Runge-Kutte est:

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 

Et je l'ai utilisé imshow() du package PyPlot pour tracer le champ u.

Answer 1

Ceci n'est pas une réponse complète, mais un avant-goût d'une tentative d'optimisation de la fonction laplacian. Le laplacian d'origine sur une matrice 10x10 m'a donné l'@Time:

0.000038 seconds (51 allocations: 12.531 KB) 

Bien que cette 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 

Donne @Time:

0.000010 seconds (6 allocations: 2.250 KB) 

Notez la réduction des allocations. Les allocations supplémentaires indiquent généralement le potentiel d'optimisation.