Hi all, very new to Julia, and not in academia BTW. I have been investigating the opinions and findings in this paper, and I thought I would contribute some Julia code which performs the aforementioned simulation to potentially ludicrous levels of accuracy, in fewer than 60 lines:
#!/usr/bin/env julia
using Printf
function main()
setprecision(trunc(Int, parse(Int, ARGS[1]) * log(10.0) / log(2.0)))
order = parse(Int, ARGS[2])
h = parse(BigFloat, ARGS[3])
steps = parse(Int, ARGS[4])
x = parse(BigFloat, ARGS[5])
y = parse(BigFloat, ARGS[6])
z = parse(BigFloat, ARGS[7])
Ļ = parse(BigFloat, ARGS[8])
Ļ = parse(BigFloat, ARGS[9])
Ī² = parse(BigFloat, ARGS[10]) / parse(BigFloat, ARGS[11])
D0 = parse(BigFloat, "0.0")
cx = [D0 for _ = 1 : order + 1]
cy = [D0 for _ = 1 : order + 1]
cz = [D0 for _ = 1 : order + 1]
wx = [D0 for _ = 1 : order]
wy = [D0 for _ = 1 : order]
wz = [D0 for _ = 1 : order]
for step = 1 : steps
@printf "%0.9e %0.9e %0.9e %0.5e\n" x y z (step - 1) * h
cx[1] = x
cy[1] = y
cz[1] = z
for k = 1 : order
wx[k] = cx[k]
wy[k] = cy[k]
wz[k] = cz[k]
cx[k+1] = Ļ * (wy[k] - wx[k]) / k
wxz = D0
for j = 1 : k
wxz += wx[j] * wz[k-j+1]
end
cy[k+1] = (Ļ * wx[k] - wxz - wy[k]) / k
wxy = D0
for j = 1 : k
wxy += wx[j] * wy[k-j+1]
end
cz[k+1] = (wxy - Ī² * wz[k]) / k
end
x = cx[order+1]
y = cy[order+1]
z = cz[order+1]
for i = order : -1 : 1
x = x * h + cx[i]
y = y * h + cy[i]
z = z * h + cz[i]
end
end
end
main()
Can anyone suggest any speedups? Would this be any use as a benchmark? BTW here is an invocation:
time -p ./tsm-lorenz.jl 160 100 .01 30001 -15.8 -17.48 35.64 10 28 8 3 >/tmp/data-Julia
[code edited to eliminate four useless lines!]
