Lorenz96をJuliaで解く その２

前回まで

ではグローバル変数を使ったJuliaのコードを書いて，実行時間0.7s程度でした。

今回はグローバル変数を使わないように，引数渡しでパラメータ設定するコードにしました。

計算時間

結果からいうと，40ms程度まで早くになりました。Fortranで計算したのと同程度の速さになりました。

コード

前回との違いは，F, N, dt, nstepをmain関数で定義して，引数で渡してるだけです。内容的には同じ計算です。

引数をそのまま書いたのでかっこわるいです。黒木さんのコードはparameter使って，かっこいいかんじになってます。https://gist.github.com/genkuroki/f0426a7f3772c804ace0e6875445ea49

function main()
    ##パラメータ
    ##Fの大きさで挙動が変わる
    F = 8
    N = 40
    dt = 0.01
    tend = 146
    nstep = Int64(tend/dt)

    xn, t =  progress(dt, F, N, nstep)
    
end

# dx/dt=f(x) 
# Lorentz96の方程式
function f(x, F, N)
    g = fill(0.0, N)
    for i=3:N-1
        g[i] = (x[i+1]-x[i-2])x[i-1] - x[i] + F
    end
    
    # 周期境界
    g[1] = (x[2]-x[N-1])x[N] - x[1] + F
    g[2] = (x[3]-x[N])x[1] - x[2] + F
    g[N] = (x[1]-x[N-2])x[N-1] - x[N] + F
    
    return g
end

# L96をRK4で解く
function RK4(xold, dt, F, N)
    k1 = f(xold, F, N)
    k2 = f(xold + k1*dt/2., F, N)
    k3 = f(xold + k2*dt/2., F, N)
    k4 = f(xold + k3*dt, F, N)
    
    xnew = xold + dt/6.0*(k1 + 2.0k2 + 2.0k3 + k4)
end

# 真値を全ステップ計算
function progress(deltat, F, N, nstep)
    xn = zeros(Float64, N, nstep)
    t = zeros(Float64, nstep)
    
    #X = fill(F,N)+randn(N)
    X = fill(Float64(F), N)
    # 初期擾乱はj=20の点にFの値の0.1%の値を与える。
    X[20] = 1.001*F
    
    for j=1:nstep
        X = RK4(X, deltat, F, N)
        for i=1:N
            xn[i, j] = X[i]
        end
        t[j] = deltat*j
    end
    
    return xn, t
end