Einführung

Ich lerne viel, weil ich die Datenassimilation für die Forschung nutzen möchte. Lorenz96-Modell, das häufig in Benchmarks zur Datenassimilation verwendet wird https://journals.ametsoc.org/jas/article/55/3/399/24564/Optimal-Sites-for-Supplementary-Weather Ich habe den Code geschrieben, um das Problem in Julia zu lösen. Ich bin mir nicht sicher, ob dies langsam oder schnell ist, also habe ich es in Python umgeschrieben und verglichen.

Modell Lorenz 96

Gleichung

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

$ J = 1 \ Punkte N $. Diesmal ist N = 40.

Wenn $ X_j = F $ anfänglich gesetzt ist, wird eine stationäre Lösung erhalten. Lösen Sie es diesmal als $ X_ {20} = 1.01F, ; F = 8 $. Wenn F = 8 ist, kann chaotisches Verhalten gesehen werden.

Dies wird mit einem 4-stufigen Präzisions-Rungekutta 4. Ordnung gelöst.

Julia Version

Die Version von Julia ist 1.5.1. Ich habe es im Jupyter-Labor ausgeführt.

using LinearAlgebra

##Parameter
##Das Verhalten ändert sich abhängig von der Größe von F.
F = 8
N = 40
dt = 0.01
tend = 146
nstep = Int32(tend/dt)

# dx/dt=f(x) 
#Lorentz 96 Gleichung
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
    
    #Periodengrenze
    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

#Löse L96 mit 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

#Berechnen Sie den wahren Wert in allen Schritten
function progress(deltat)
    xn = zeros(Float64, N, nstep)
    t = zeros(Float64, nstep)
    
    #X = fill(F,N)+randn(N)
    X = fill(Float64(F), N)
    #Die anfängliche Störung ist j=0 von F-Wert bei 20 Punkten.1%Gibt den Wert von.
    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_Berechnung von wahr
@time xn, t = progress(dt)

Python-Version

Python ist 3.7.3. Python berechnet mit numpy, das mit pip installiert ist. Numpy in Anaconda ist schneller, aber es ist vorerst ein Test, also ist es mir egal.

import numpy as np
from copy import copy

##Parameter
##Das Verhalten ändert sich abhängig von der Größe von F.
F = 8
#Anzahl der Analysedatenpunkte (Dimension)
N = 40
dt = 0.01
#1 Jahr Berechnung
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

#Löse L96 mit 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

#Berechnen Sie den wahren Wert in allen Schritten
def progress(deltat):
    xn = np.zeros((N, nstep))
    t = np.zeros((nstep))
    
    X = float(F)*np.ones(N)
    #X = fill(Float64(F), N)
    #Die anfängliche Störung ist j=0 des Wertes von F am Punkt 19 (unter Berücksichtigung des Beginns von 0).1%Gibt den Wert von.
    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

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

Berechnungszeit

Julia 0.73 s Python 2.94 s

Und Julia war schneller. 2,3 GHz Quad-Core Intel Core i5 Ich habe es mit einem geladenen MacBook Pro gemacht.

Ich habe nicht wirklich darüber nachgedacht, Python zu optimieren, daher denke ich, dass es schneller sein wird. Ich habe die Reihenfolge der Array-Indizes geändert, aber es wurde nicht schneller. Ich kann JIT-Kompilierung mit Python durchführen, aber ich habe es nicht getan, also habe ich es dieses Mal auf Eis gelegt.

Vorerst schien Julias Rechengeschwindigkeit nicht langsam zu sein. Übrigens habe ich auch den Code für die Datenassimilation durch Extended Kalman Filter geschrieben, aber Julia war schneller.

Bonus: Zeichnen

Code zum Zeichnen des Berechnungsergebnisses in Julia

Hoffmeler-Diagramm

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")

Zeitumstellung

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: Datenspeicherung

Als HDF5 speichern. Speichern Sie die Daten alle 5 Schritte in der zweiten Hälfte des Berechnungsergebnisses.

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 sieht so aus

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

Lösen des Lorenz 96-Modells mit Julia und Python