introduction

J'étudie beaucoup parce que je veux utiliser l'assimilation de données pour la recherche. Modèle Lorenz96 souvent utilisé dans les benchmarks d'assimilation de données https://journals.ametsoc.org/jas/article/55/3/399/24564/Optimal-Sites-for-Supplementary-Weather J'ai écrit le code pour résoudre le problème de Julia. Je ne sais pas si c'est lent ou rapide, alors je l'ai réécrit en Python et je l'ai comparé.

Modèle Lorenz 96

équation

\frac{dX_j}{dt} = (X_{j+1} - X_{j-2})X_{j-1} - X_j +F

Cependant, $ j = 1 \ points N $. Cette fois N = 40.

Si $ X_j = F $ est initialisé, une solution stationnaire sera obtenue. Cette fois, résolvez-le comme $ X_ {20} = 1.01F, ; F = 8 $. Lorsque F = 8, un comportement chaotique peut être observé.

Ce problème est résolu avec un Rungekutta de précision de 4ème ordre en 4 étapes.

Version Julia

La version de Julia est 1.5.1. Je l'ai exécuté dans le laboratoire Jupyter.

using LinearAlgebra

##Paramètres
##Changements de comportement en fonction de la taille de F
F = 8
N = 40
dt = 0.01
tend = 146
nstep = Int32(tend/dt)

# dx/dt=f(x) 
#Équation de Lorentz 96
function f(x)
    f = fill(0.0, N)
    for i=3:N-1
        f[i] = (x[i+1]-x[i-2])x[i-1] - x[i] + F
    end
    
    #Limite de période
    f[1] = (x[2]-x[N-1])x[N] - x[1] + F
    f[2] = (x[3]-x[N])x[1] - x[2] + F
    f[N] = (x[1]-x[N-2])x[N-1] - x[N] + F
    
    return f
end

#Résoudre L96 avec RK4
function RK4(xold, dt)
    k1 = f(xold)
    k2 = f(xold + k1*dt/2.)
    k3 = f(xold + k2*dt/2.)
    k4 = f(xold + k3*dt)
    
    xnew = xold + dt/6.0*(k1 + 2.0k2 + 2.0k3 + k4)
end

#Calculez la vraie valeur dans toutes les étapes
function progress(deltat)
    xn = zeros(Float64, N, nstep)
    t = zeros(Float64, nstep)
    
    #X = fill(F,N)+randn(N)
    X = fill(Float64(F), N)
    #La perturbation initiale est j=0 de la valeur F à 20 points.1%Donne la valeur de.
    X[20] = 1.001*F
    
    for j=1:nstep
        X = RK4(X, deltat)
        for i=1:N
            xn[i, j] = X[i]
        end
        t[j] = dt*j
    end
    
    return xn, t
end

# x_Calcul de vrai
@time xn, t = progress(dt)

Version Python

Python est 3.7.3. Python calcule avec numpy installé avec pip. Numpy dans Anaconda est plus rapide, mais c'est un test pour le moment, donc je m'en fiche.

import numpy as np
from copy import copy

##Paramètres
##Changements de comportement en fonction de la taille de F
F = 8
#Nombre de points de données d'analyse (dimension)
N = 40
dt = 0.01
#Calcul sur 1 an
tend=146
nstep = int(tend/dt)

def f(x):
    f = np.zeros((N))
    for i in range(2, N-1):
        f[i] = (x[i+1]-x[i-2])*x[i-1] - x[i] + F
    
    f[0] = (x[1]-x[N-2])*x[N-1] - x[0] + F
    f[1] = (x[2]-x[N-1])*x[0] - x[1] + F
    f[N-1] = (x[0]-x[N-3])*x[N-2] - x[N-1] + F
    
    return f

#Résoudre L96 avec RK4
def Model(xold, dt):
    k1 = f(xold)
    k2 = f(xold + k1*dt/2.)
    k3 = f(xold + k2*dt/2.)
    k4 = f(xold + k3*dt)
    
    xnew = xold + dt/6.0*(k1 + 2.0*k2 + 2.0*k3 + k4)
    return xnew

#Calculez la vraie valeur dans toutes les étapes
def progress(deltat):
    xn = np.zeros((N, nstep))
    t = np.zeros((nstep))
    
    X = float(F)*np.ones(N)
    #X = fill(Float64(F), N)
    #La perturbation initiale est j=0 de la valeur de F au point de 19 (en considérant le début de 0).1%Donne la valeur de.
    X[19] = 1.001*F
    
    for i in range(N):
        xn[i, 0] = copy(X[i])
    
    for j in range(1, nstep):
        X = Model(X, deltat)
        xn[:, j] = X[:]
        #print(X[20], X[1], deltat)
        t[j] = dt*j

    return xn, t

#Courir
start = time.time()
xn, t = progress(dt)
elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

Temps de calcul

Julia 0.73 s Python 2.94 s

Et Julia était plus rapide. Intel Core i5 quadricœur 2,3 GHz Je l'ai fait avec un MacBook Pro chargé.

Je n'ai pas vraiment pensé à optimiser Python, donc je pense que ce sera plus rapide. J'ai changé l'ordre des index du tableau, mais cela ne s'est pas accéléré. Je peux faire une compilation JIT avec Python, mais je ne l'ai pas fait, donc je l'ai mise en attente cette fois.

Pour le moment, la vitesse de calcul de Julia ne semblait pas lente. Au fait, j'ai également écrit le code pour l'assimilation des données par Extended Kalman Filter, mais Julia était plus rapide.

Bonus: dessin

Code pour tracer le résultat du calcul dans Julia

Diagramme de Hoffmeler

using PyPlot

fig = plt.figure(figsize=(6, 5.))
ax1 = fig.add_subplot(111)

ax1.set_title("true")
ax1.set_xlabel("X_i")
ax1.set_ylabel("time step")

image1 = ax1.contourf(xn'[1:5000,:], )
cb1 = fig.colorbar(image1, ax=ax1, label="X")

changement d'heure

using PyPlot

fig = plt.figure(figsize=(5, 8.))
ax1 = fig.add_subplot(211)
nhalf = Int(nstep/2)
nstart = nhalf
nend = nstep
ax1.plot(t[nstart:nend], xn[1, nstart:nend], label="x1")
ax1.plot(t[nstart:nend], xn[2, nstart:nend], label="x2")
ax1.plot(t[nstart:nend], xn[3, nstart:nend], label="x3")

Bonus: stockage de données

Enregistrez sous HDF5. Sauvegardez les données toutes les 5 étapes dans la seconde moitié du résultat du calcul.

using HDF5

h5open("L96_true_obs.h5", "w") do file
    write(file, "Xn_true", xn[:, nhalf:5:nstep])
    write(file, "t", t[nhalf:5:nstep])
end

Reed ressemble à ça

file = h5open("L96_true_obs.h5", "r") 
Xn_true = read(file, "Xn_true")
t = read(file, "t")
close(file)

Résolution du modèle Lorenz 96 avec Julia et Python