Merged
Conversation
- Extend DynamicalModels with parameter management via ParamIO + DataVault: configs/ (lorenz, rossler, logistic with inheritance) demonstrate 3 paramset patterns (plain path_keys, multi-block union, config inheritance) - Add docs/src/tutorials/paramset.jl: Literate script executed at doc build time, embeds Lyapunov exponent sweeps and bifurcation diagram plots as inline images - Update docs/make.jl: add favicon/logo from GitHub profile, MathJax3, canonical URL, Literate execution loop, proper deploydocs settings - Fix LogisticMap functor signature to (n, x) to match map_solver interface - Fix docs/Project.toml: use relative path ".." for DynamicalModels source - Update Documentation.yml: clone DataVault.jl and ParamIO.jl as sibling repos so their local [sources] paths resolve correctly in CI - Add integration test for data/code separation (test/integration/test_separation.jl) - Update .gitignore: exclude docs/build/, docs/src/assets/, generated Literate .md, and untrack docs/build/ and root Manifest.toml Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
- [sources] in Project.toml and docs/Project.toml now point to GitHub URLs instead of local sibling paths so CI can resolve them without cloning - Remove clone-sibling-repos step from Documentation.yml - Add docs/Manifest.toml to .gitignore; local env resolves via path but CI generates fresh from GitHub URLs - Bump version to v0.2.0 Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Contributor
|
Your PR requires formatting changes to meet the project's style guidelines. Click here to view the suggested changes.diff --git a/Project.toml b/Project.toml
index ddb6fc4..3ba9283 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,16 +15,18 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-[sources]
-DataVault = {url = "https://github.com/sotashimozono/DataVault.jl.git"}
-ParamIO = {url = "https://github.com/sotashimozono/ParamIO.jl.git"}
+[sources.DataVault]
+url = "https://github.com/sotashimozono/DataVault.jl.git"
+
+[sources.ParamIO]
+url = "https://github.com/sotashimozono/ParamIO.jl.git"
[compat]
-FFTW = "1.10.0"
-ForwardDiff = "1.3.0"
-LinearAlgebra = "1.12.0"
-Literate = "2.21.0"
+FFTW = "1.10"
+ForwardDiff = "1.3"
+LinearAlgebra = "1.12"
+Literate = "2.21"
Plots = "1.41.2"
-Random = "1.11.0"
+Random = "1.11"
Statistics = "1.11.1"
-Test = "1.11.0"
+Test = "1.11"
diff --git a/docs/Project.toml b/docs/Project.toml
index 6ef5157..4ec8c96 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -7,7 +7,11 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
ParamIO = "938a3ac2-d340-473c-bcf1-88af577e4ccf"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-[sources]
-DataVault = {path = "../../DataVault.jl"}
-DynamicalModels = {path = "/home/souta/work/dev/DynamicalModels.jl"}
-ParamIO = {path = "../../ParamIO.jl"}
+[sources.DataVault]
+path = "../../DataVault.jl"
+
+[sources.DynamicalModels]
+path = "/home/souta/work/dev/DynamicalModels.jl"
+
+[sources.ParamIO]
+path = "../../ParamIO.jl"
diff --git a/docs/make.jl b/docs/make.jl
index c886d67..c8bebe0 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,13 +1,16 @@
-using DynamicalModels
using Documenter
-using Literate
using Downloads
+using DynamicalModels
+using Literate
# ── Assets: favicon / logo from GitHub profile ────────────────────────────────
assets_dir = joinpath(@__DIR__, "src", "assets")
mkpath(assets_dir)
-Downloads.download("https://github.com/sotashimozono.png", joinpath(assets_dir, "favicon.ico"))
+Downloads.download(
+ "https://github.com/sotashimozono.png",
+ joinpath(assets_dir, "favicon.ico")
+)
Downloads.download("https://github.com/sotashimozono.png", joinpath(assets_dir, "logo.png"))
# ── Literate: tutorial scripts → markdown ─────────────────────────────────────
@@ -20,48 +23,48 @@ for script in filter(f -> endswith(f, ".jl"), readdir(TUTORIALS_SRC))
joinpath(TUTORIALS_SRC, script),
TUTORIALS_OUT;
documenter = true,
- execute = true,
- credit = false,
+ execute = true,
+ credit = false
)
end
# ── Documenter ────────────────────────────────────────────────────────────────
-makedocs(
+makedocs(;
sitename = "DynamicalModels.jl",
- modules = [DynamicalModels],
- format = Documenter.HTML(
- canonical = "https://codes.sota-shimozono.com/DynamicalModels.jl/stable/",
+ modules = [DynamicalModels],
+ format = Documenter.HTML(;
+ canonical = "https://codes.sota-shimozono.com/DynamicalModels.jl/stable/",
prettyurls = get(ENV, "CI", "false") == "true",
mathengine = MathJax3(
Dict(
:tex => Dict(
:inlineMath => [["\$", "\$"], ["\\(", "\\)"]],
- :tags => "ams",
- :packages => ["base", "ams", "autoload", "physics"],
- ),
- ),
+ :tags => "ams",
+ :packages => ["base", "ams", "autoload", "physics"]
+ )
+ )
),
- assets = ["assets/favicon.ico"],
+ assets = ["assets/favicon.ico"]
),
pages = [
"Home" => "index.md",
"Models" => [
"Van der Pol Oscillator" => "models/vanderpol.md",
- "Lorenz System" => "models/lorenz.md",
- "Rössler System" => "models/rossler.md",
- "Hodgkin-Huxley Model" => "models/hodgkin-huxley.md",
+ "Lorenz System" => "models/lorenz.md",
+ "Rössler System" => "models/rossler.md",
+ "Hodgkin-Huxley Model" => "models/hodgkin-huxley.md",
],
"Analysis Tools" => "analysis.md",
"Tutorials" => [
"Parameter Space Management" => "tutorials/paramset.md",
],
- ],
+ ]
)
-deploydocs(
- repo = "github.com/sotashimozono/DynamicalModels.jl.git",
- devbranch = "main",
- target = "build",
- push_preview = false,
+deploydocs(;
+ repo = "github.com/sotashimozono/DynamicalModels.jl.git",
+ devbranch = "main",
+ target = "build",
+ push_preview = false
)
diff --git a/docs/src/tutorials/paramset.jl b/docs/src/tutorials/paramset.jl
index 0118f6e..8e9620a 100644
--- a/docs/src/tutorials/paramset.jl
+++ b/docs/src/tutorials/paramset.jl
@@ -9,10 +9,14 @@
ENV["GKSwstype"] = "100" # headless Plots rendering
-using DynamicalModels, DataVault, ParamIO, Printf, Plots
+using DataVault
+using DynamicalModels
+using ParamIO
+using Plots
+using Printf
const CONFIGS = joinpath(@__DIR__, "..", "..", "..", "configs")
-const OUTDIR = mktempdir() # data lives here; thrown away after docs build
+const OUTDIR = mktempdir() # data lives here; thrown away after docs build
# ---
# ## Pattern 1 — Plain `path_keys`: Lorenz ρ sweep
@@ -28,28 +32,36 @@ spec_lorenz = ParamIO.load(joinpath(CONFIGS, "lorenz.toml"))
println("path_keys : ", spec_lorenz.path_keys)
println("DataKeys : ", length(ParamIO.expand(spec_lorenz)), " (param points × samples)")
-vault_lorenz = Vault(joinpath(CONFIGS, "lorenz.toml"); outdir=joinpath(OUTDIR, "lorenz"))
+vault_lorenz = Vault(joinpath(CONFIGS, "lorenz.toml"); outdir = joinpath(OUTDIR, "lorenz"))
-for key in DataVault.keys(vault_lorenz; status=:pending)
- ρ = Float64(key.params["system.rho"])
- model = Lorenz(; σ=10.0, ρ=ρ, β=8/3)
- x0 = [1.0, 1.0, 1.0]
- λ = lyapunov_exponent(model, x0, 0.1; warmup=300, n_iterations=2000)
+for key in DataVault.keys(vault_lorenz; status = :pending)
+ ρ = Float64(key.params["system.rho"])
+ model = Lorenz(; σ = 10.0, ρ = ρ, β = 8 / 3)
+ x0 = [1.0, 1.0, 1.0]
+ λ = lyapunov_exponent(model, x0, 0.1; warmup = 300, n_iterations = 2000)
DataVault.save!(vault_lorenz, key, Dict("rho" => ρ, "lyapunov" => λ))
- mark_done!(vault_lorenz, key; tag_value=λ)
+ mark_done!(vault_lorenz, key; tag_value = λ)
end
println("\n| ρ | λ (Lyapunov) | chaotic? |")
println("|-------|--------------|----------|")
seen = Set{Float64}()
-rho_vals = Float64[]
-lam_vals = Float64[]
-for key in sort(DataVault.keys(vault_lorenz; status=:done); by=k->k.params["system.rho"])
+rho_vals = Float64[]
+lam_vals = Float64[]
+for key in
+ sort(DataVault.keys(vault_lorenz; status = :done); by = k -> k.params["system.rho"])
ρ = Float64(key.params["system.rho"])
- ρ ∈ seen && continue; push!(seen, ρ)
+ ρ ∈ seen && continue
+ push!(seen, ρ)
d = DataVault.load(vault_lorenz, key)
- push!(rho_vals, d["rho"]); push!(lam_vals, d["lyapunov"])
- @printf("| %5.2f | %12.4f | %-8s |\n", d["rho"], d["lyapunov"], d["lyapunov"] > 0 ? "yes" : "no")
+ push!(rho_vals, d["rho"])
+ push!(lam_vals, d["lyapunov"])
+ @printf(
+ "| %5.2f | %12.4f | %-8s |\n",
+ d["rho"],
+ d["lyapunov"],
+ d["lyapunov"] > 0 ? "yes" : "no"
+ )
end
# The output path for each key:
@@ -62,13 +74,14 @@ println("(plain name → no group prefix)")
# The largest Lyapunov exponent crosses zero at ρ ≈ 24.74, marking the onset
# of chaos. Positive λ → exponential divergence of nearby trajectories.
-p_lorenz = scatter(rho_vals, lam_vals;
+p_lorenz = scatter(
+ rho_vals, lam_vals;
xlabel = "ρ", ylabel = "λ (Lyapunov exponent)",
- title = "Lorenz system — λ vs ρ",
+ title = "Lorenz system — λ vs ρ",
legend = false,
- marker = (:circle, 8),
+ marker = (:circle, 8)
)
-hline!(p_lorenz, [0.0]; linestyle=:dash, color=:gray, label="λ = 0")
+hline!(p_lorenz, [0.0]; linestyle = :dash, color = :gray, label = "λ = 0")
annotate!(p_lorenz, 24.74, 0.05, text("ρ ≈ 24.74", 9, :left))
p_lorenz
@@ -84,49 +97,71 @@ p_lorenz
spec_rossler = ParamIO.load(joinpath(CONFIGS, "rossler.toml"))
println("\nRössler union: $(length(spec_rossler.paramsets)) blocks")
-println("Total points : $(length(ParamIO.expand(spec_rossler)) ÷ spec_rossler.study.total_samples) param points")
-println("(vs Cartesian: $(5*5) points — most would be off the bifurcation paths)")
+println(
+ "Total points : $(length(ParamIO.expand(spec_rossler)) ÷ spec_rossler.study.total_samples) param points"
+)
+println("(vs Cartesian: $(5 * 5) points — most would be off the bifurcation paths)")
-vault_rossler = Vault(joinpath(CONFIGS, "rossler.toml"); outdir=joinpath(OUTDIR, "rossler"))
+vault_rossler =
+ Vault(joinpath(CONFIGS, "rossler.toml"); outdir = joinpath(OUTDIR, "rossler"))
-for key in DataVault.keys(vault_rossler; status=:pending)
+for key in DataVault.keys(vault_rossler; status = :pending)
a = Float64(key.params["system.a"])
b = Float64(key.params["system.b"])
c = Float64(key.params["system.c"])
- model = Rossler(; a=a, b=b, c=c)
- x0 = [1.0, 1.0, 1.0]
- λ = lyapunov_exponent(model, x0, 0.1; warmup=200, n_iterations=1500)
+ model = Rossler(; a = a, b = b, c = c)
+ x0 = [1.0, 1.0, 1.0]
+ λ = lyapunov_exponent(model, x0, 0.1; warmup = 200, n_iterations = 1500)
DataVault.save!(vault_rossler, key, Dict("a" => a, "b" => b, "c" => c, "lyapunov" => λ))
- mark_done!(vault_rossler, key; tag_value=λ)
+ mark_done!(vault_rossler, key; tag_value = λ)
end
-a_vals_r1 = Float64[]; lam_vals_r1 = Float64[]
-c_vals_r2 = Float64[]; lam_vals_r2 = Float64[]
+a_vals_r1 = Float64[];
+lam_vals_r1 = Float64[]
+c_vals_r2 = Float64[];
+lam_vals_r2 = Float64[]
println("\nBlock 1 (c=14.0, a sweep — period-doubling route):")
println("| a | λ | chaotic? |")
println("|------|--------|----------|")
seen_r = Set{Tuple}()
-for key in sort(DataVault.keys(vault_rossler; status=:done); by=k->(k.params["system.c"],k.params["system.a"]))
+for key in sort(
+ DataVault.keys(vault_rossler; status = :done);
+ by = k -> (k.params["system.c"], k.params["system.a"])
+ )
d = DataVault.load(vault_rossler, key)
- t = (round(d["a"],digits=3), round(d["c"],digits=3))
- t ∈ seen_r && continue; push!(seen_r, t)
+ t = (round(d["a"]; digits = 3), round(d["c"]; digits = 3))
+ t ∈ seen_r && continue
+ push!(seen_r, t)
abs(d["c"] - 14.0) < 0.01 || continue
- push!(a_vals_r1, d["a"]); push!(lam_vals_r1, d["lyapunov"])
- @printf("| %.2f | %6.3f | %-8s |\n", d["a"], d["lyapunov"], d["lyapunov"] > 0 ? "yes" : "no")
+ push!(a_vals_r1, d["a"])
+ push!(lam_vals_r1, d["lyapunov"])
+ @printf(
+ "| %.2f | %6.3f | %-8s |\n",
+ d["a"],
+ d["lyapunov"],
+ d["lyapunov"] > 0 ? "yes" : "no"
+ )
end
println("\nBlock 2 (a=0.2, c sweep — spiral→screw transition):")
println("| c | λ | chaotic? |")
println("|------|--------|----------|")
seen_r2 = Set{Tuple}()
-for key in sort(DataVault.keys(vault_rossler; status=:done); by=k->k.params["system.c"])
+for key in sort(DataVault.keys(vault_rossler; status = :done); by = k -> k.params["system.c"])
d = DataVault.load(vault_rossler, key)
- t = (round(d["a"],digits=3), round(d["c"],digits=3))
- t ∈ seen_r2 && continue; push!(seen_r2, t)
+ t = (round(d["a"]; digits = 3), round(d["c"]; digits = 3))
+ t ∈ seen_r2 && continue
+ push!(seen_r2, t)
abs(d["a"] - 0.2) < 0.01 && abs(d["c"] - 14.0) > 0.1 || continue
- push!(c_vals_r2, d["c"]); push!(lam_vals_r2, d["lyapunov"])
- @printf("| %.1f | %6.3f | %-8s |\n", d["c"], d["lyapunov"], d["lyapunov"] > 0 ? "yes" : "no")
+ push!(c_vals_r2, d["c"])
+ push!(lam_vals_r2, d["lyapunov"])
+ @printf(
+ "| %.1f | %6.3f | %-8s |\n",
+ d["c"],
+ d["lyapunov"],
+ d["lyapunov"] > 0 ? "yes" : "no"
+ )
end
# ### Rössler: Lyapunov exponents along two bifurcation routes
@@ -135,21 +170,23 @@ end
# A Cartesian product would mix these routes, producing 25 physically irrelevant
# combinations. The union design captures exactly the two bifurcation paths.
-p_r1 = scatter(a_vals_r1, lam_vals_r1;
+p_r1 = scatter(
+ a_vals_r1, lam_vals_r1;
xlabel = "a (c = 14.0 fixed)", ylabel = "λ",
- title = "Rössler — period-doubling route",
- legend = false, marker = (:circle, 8),
+ title = "Rössler — period-doubling route",
+ legend = false, marker = (:circle, 8)
)
-hline!(p_r1, [0.0]; linestyle=:dash, color=:gray)
+hline!(p_r1, [0.0]; linestyle = :dash, color = :gray)
-p_r2 = scatter(c_vals_r2, lam_vals_r2;
+p_r2 = scatter(
+ c_vals_r2, lam_vals_r2;
xlabel = "c (a = 0.2 fixed)", ylabel = "λ",
- title = "Rössler — spiral→screw transition",
- legend = false, marker = (:diamond, 8),
+ title = "Rössler — spiral→screw transition",
+ legend = false, marker = (:diamond, 8)
)
-hline!(p_r2, [0.0]; linestyle=:dash, color=:gray)
+hline!(p_r2, [0.0]; linestyle = :dash, color = :gray)
-plot(p_r1, p_r2; layout=(1,2), size=(820, 360))
+plot(p_r1, p_r2; layout = (1, 2), size = (820, 360))
# ---
# ## Pattern 3 — Config inheritance: Logistic map bifurcation diagram
@@ -158,43 +195,56 @@ plot(p_r1, p_r2; layout=(1,2), size=(820, 360))
# This lets you run the base study first, then extend without recomputing.
spec_default = ParamIO.load(joinpath(CONFIGS, "logistic", "default.toml"))
-spec_chaos = ParamIO.load(joinpath(CONFIGS, "logistic", "chaos.toml"))
-println("\nLogistic default : $(length(ParamIO.expand(spec_default))) DataKeys ($(spec_default.study.total_samples) sample)")
-println("Logistic chaos : $(length(ParamIO.expand(spec_chaos))) DataKeys (inherits default + adds chaos regime)")
+spec_chaos = ParamIO.load(joinpath(CONFIGS, "logistic", "chaos.toml"))
+println(
+ "\nLogistic default : $(length(ParamIO.expand(spec_default))) DataKeys ($(spec_default.study.total_samples) sample)"
+)
+println(
+ "Logistic chaos : $(length(ParamIO.expand(spec_chaos))) DataKeys (inherits default + adds chaos regime)"
+)
-vault_logistic = Vault(joinpath(CONFIGS, "logistic", "chaos.toml"); outdir=joinpath(OUTDIR, "logistic"))
+vault_logistic = Vault(
+ joinpath(CONFIGS, "logistic", "chaos.toml");
+ outdir = joinpath(OUTDIR, "logistic")
+)
-for key in DataVault.keys(vault_logistic; status=:pending)
- r = Float64(key.params["system.r"])
+for key in DataVault.keys(vault_logistic; status = :pending)
+ r = Float64(key.params["system.r"])
n_iter = Int(key.params["system.n_iter"])
n_skip = Int(key.params["system.n_skip"])
- model = LogisticMap(; r=r)
- x0 = [0.5]
- raw = map_solver(model, n_iter, x0) # (n_iter × 1) matrix
- orbit = raw[:, 1] # 1-D time series
- attractor = orbit[n_skip+1:end]
+ model = LogisticMap(; r = r)
+ x0 = [0.5]
+ raw = map_solver(model, n_iter, x0) # (n_iter × 1) matrix
+ orbit = raw[:, 1] # 1-D time series
+ attractor = orbit[(n_skip + 1):end]
- uniq = unique(round.(attractor; digits=4))
+ uniq = unique(round.(attractor; digits = 4))
period = length(uniq) > 50 ? Inf : length(uniq)
- DataVault.save!(vault_logistic, key, Dict(
- "r" => r,
- "period" => period,
- "attractor" => attractor,
- ))
- mark_done!(vault_logistic, key; tag_value=Float64(period))
+ DataVault.save!(
+ vault_logistic,
+ key,
+ Dict(
+ "r" => r,
+ "period" => period,
+ "attractor" => attractor
+ )
+ )
+ mark_done!(vault_logistic, key; tag_value = Float64(period))
end
println("\n| r | period (approx) | regime |")
println("|------|-----------------|---------------|")
seen_l = Set{Float64}()
-for key in sort(DataVault.keys(vault_logistic; status=:done); by=k->k.params["system.r"])
+for key in
+ sort(DataVault.keys(vault_logistic; status = :done); by = k -> k.params["system.r"])
r = Float64(key.params["system.r"])
- r ∈ seen_l && continue; push!(seen_l, r)
- d = DataVault.load(vault_logistic, key)
- p = d["period"]
- reg = isinf(p) ? "chaotic" : p == 1 ? "fixed point" : "period-$(Int(p))"
+ r ∈ seen_l && continue
+ push!(seen_l, r)
+ d = DataVault.load(vault_logistic, key)
+ p = d["period"]
+ reg = isinf(p) ? "chaotic" : p == 1 ? "fixed point" : "period-$(Int(p))"
@printf("| %.2f | %-15s | %-13s |\n", r, isinf(p) ? "∞" : string(Int(p)), reg)
end
@@ -203,26 +253,29 @@ end
# Each vertical column of dots shows the long-run attractor for a given `r`.
# The period-doubling cascade and the onset of chaos at `r ≈ 3.57` are clearly visible.
-r_plot = Float64[]
+r_plot = Float64[]
att_plot = Float64[]
-for key in sort(DataVault.keys(vault_logistic; status=:done); by=k->k.params["system.r"])
+for key in
+ sort(DataVault.keys(vault_logistic; status = :done); by = k -> k.params["system.r"])
r = Float64(key.params["system.r"])
r ∈ Set(r_plot) && continue
d = DataVault.load(vault_logistic, key)
for v in d["attractor"]
- push!(r_plot, r); push!(att_plot, v)
+ push!(r_plot, r)
+ push!(att_plot, v)
end
end
-scatter(r_plot, att_plot;
- xlabel = "r",
- ylabel = "x (attractor)",
- title = "Logistic map — bifurcation diagram",
+scatter(
+ r_plot, att_plot;
+ xlabel = "r",
+ ylabel = "x (attractor)",
+ title = "Logistic map — bifurcation diagram",
markersize = 1,
markerstrokewidth = 0,
- legend = false,
- alpha = 0.3,
- size = (700, 400),
+ legend = false,
+ alpha = 0.3,
+ size = (700, 400)
)
# ---
diff --git a/examples/lorenz_analysis.jl b/examples/lorenz_analysis.jl
index 69d8a8c..15b16ef 100644
--- a/examples/lorenz_analysis.jl
+++ b/examples/lorenz_analysis.jl
@@ -15,7 +15,7 @@ println("Lorenz System Analysis")
println("="^60)
# Create Lorenz model with standard parameters
-model = Lorenz(σ=10.0, ρ=28.0, β=8/3)
+model = Lorenz(; σ = 10.0, ρ = 28.0, β = 8 / 3)
println("\nModel parameters:")
println(" σ = $(model.σ)")
println(" ρ = $(model.ρ)")
@@ -31,12 +31,12 @@ println("-"^60)
# Calculate largest Lyapunov exponent
println("\nCalculating largest Lyapunov exponent...")
-λ_max = lyapunov_exponent(model, x0, 0.1; warmup=1000, n_iterations=5000)
+λ_max = lyapunov_exponent(model, x0, 0.1; warmup = 1000, n_iterations = 5000)
println("Largest Lyapunov exponent: ", λ_max)
if λ_max > 0
println("✓ System is CHAOTIC (λ > 0)")
- println(" Lyapunov time (predictability horizon): ", 1/λ_max, " time units")
+ println(" Lyapunov time (predictability horizon): ", 1 / λ_max, " time units")
else
println("✗ System is NOT chaotic")
end
@@ -44,7 +44,7 @@ end
# Calculate full Lyapunov spectrum
println("\nCalculating full Lyapunov spectrum...")
println("(This may take a while...)")
-λs = lyapunov_spectrum(model, x0, 0.5; warmup=1000, n_iterations=2000)
+λs = lyapunov_spectrum(model, x0, 0.5; warmup = 1000, n_iterations = 2000)
println("\nLyapunov spectrum:")
for (i, λ) in enumerate(λs)
println(" λ[$i] = ", λ)
@@ -65,14 +65,18 @@ plane_normal = [0.0, 0.0, 1.0]
plane_point = [0.0, 0.0, 27.0]
println("\nCalculating Poincaré section at z = 27...")
-section = poincare_section(model, x0, plane_normal, plane_point, 500.0;
- dt=0.01, direction=:both)
+section = poincare_section(
+ model, x0, plane_normal, plane_point, 500.0;
+ dt = 0.01, direction = :both
+)
println("Number of section crossings: ", length(section))
# Get 2D projection
-x_coords, y_coords = poincare_map_2d(model, x0, (1, 2),
- plane_normal, plane_point, 500.0)
+x_coords, y_coords = poincare_map_2d(
+ model, x0, (1, 2),
+ plane_normal, plane_point, 500.0
+)
println("\nPoincaré section statistics:")
println(" x range: [", minimum(x_coords), ", ", maximum(x_coords), "]")
@@ -112,8 +116,8 @@ println("Analysis Complete")
println("="^60)
println("\nSummary:")
println(" - Lorenz system exhibits chaotic behavior")
-println(" - Positive Lyapunov exponent: ", round(λ_max, digits=3))
-println(" - Fractal dimension: ", round(D_KY, digits=3))
+println(" - Positive Lyapunov exponent: ", round(λ_max; digits = 3))
+println(" - Fractal dimension: ", round(D_KY; digits = 3))
println(" - Poincaré section reveals strange attractor structure")
println(" - System is dissipative (volume contracts)")
println("\n" * "="^60)
diff --git a/examples/rossler_analysis.jl b/examples/rossler_analysis.jl
index edcf75e..efea8d5 100644
--- a/examples/rossler_analysis.jl
+++ b/examples/rossler_analysis.jl
@@ -14,7 +14,7 @@ println("Rössler System Analysis")
println("="^60)
# Create Rössler model with standard chaotic parameters
-model = Rossler(a=0.2, b=0.2, c=5.7)
+model = Rossler(; a = 0.2, b = 0.2, c = 5.7)
println("\nModel parameters:")
println(" a = $(model.a)")
println(" b = $(model.b)")
@@ -30,7 +30,7 @@ println("-"^60)
# Calculate largest Lyapunov exponent
println("\nCalculating largest Lyapunov exponent...")
-λ_max = lyapunov_exponent(model, x0, 0.1; warmup=1000, n_iterations=5000)
+λ_max = lyapunov_exponent(model, x0, 0.1; warmup = 1000, n_iterations = 5000)
println("Largest Lyapunov exponent: ", λ_max)
if λ_max > 0
@@ -50,15 +50,19 @@ plane_normal = [0.0, 1.0, 0.0]
plane_point = [0.0, 0.0, 0.0]
println("\nCalculating Poincaré section at y = 0...")
-section = poincare_section(model, x0, plane_normal, plane_point, 500.0;
- dt=0.01, direction=:positive)
+section = poincare_section(
+ model, x0, plane_normal, plane_point, 500.0;
+ dt = 0.01, direction = :positive
+)
println("Number of section crossings: ", length(section))
# Get x-z projection
-x_coords, z_coords = poincare_map_2d(model, x0, (1, 3),
- plane_normal, plane_point, 500.0;
- direction=:positive)
+x_coords, z_coords = poincare_map_2d(
+ model, x0, (1, 3),
+ plane_normal, plane_point, 500.0;
+ direction = :positive
+)
println("\nPoincaré section statistics:")
println(" x range: [", minimum(x_coords), ", ", maximum(x_coords), "]")
@@ -67,9 +71,9 @@ println(" z range: [", minimum(z_coords), ", ", maximum(z_coords), "]")
# Return map analysis
if length(x_coords) > 1
println("\nReturn map analysis (x-coordinates):")
- x_n = x_coords[1:end-1]
+ x_n = x_coords[1:(end - 1)]
x_n1 = x_coords[2:end]
-
+
println(" Number of points in return map: ", length(x_n))
println(" Return map reveals structure of attractor")
println(" (Plot x_n vs x_n1 to see return map)")
@@ -86,9 +90,9 @@ println()
c_values = [2.0, 3.0, 4.0, 4.5, 5.0, 5.5, 5.7, 6.0]
for c in c_values
- test_model = Rossler(a=0.2, b=0.2, c=c)
- λ = lyapunov_exponent(test_model, x0, 0.1; warmup=500, n_iterations=1000)
-
+ test_model = Rossler(; a = 0.2, b = 0.2, c = c)
+ λ = lyapunov_exponent(test_model, x0, 0.1; warmup = 500, n_iterations = 1000)
+
behavior = if λ > 0.01
"CHAOTIC"
elseif abs(λ) < 0.01
@@ -96,8 +100,8 @@ for c in c_values
else
"STABLE"
end
-
- println(" c = ", c, "\tλ = ", round(λ, digits=4), "\t[", behavior, "]")
+
+ println(" c = ", c, "\tλ = ", round(λ; digits = 4), "\t[", behavior, "]")
end
println("\nObservations:")
@@ -128,7 +132,7 @@ println("Analysis Complete")
println("="^60)
println("\nSummary:")
println(" - Rössler system exhibits chaotic behavior at c = 5.7")
-println(" - Lyapunov exponent: ", round(λ_max, digits=3))
+println(" - Lyapunov exponent: ", round(λ_max; digits = 3))
println(" - Poincaré section shows characteristic bun shape")
println(" - Period-doubling route to chaos observed")
println(" - Ideal for studying transition to chaos")
diff --git a/examples/vanderpol_analysis.jl b/examples/vanderpol_analysis.jl
index f800e4f..d88f460 100644
--- a/examples/vanderpol_analysis.jl
+++ b/examples/vanderpol_analysis.jl
@@ -18,7 +18,7 @@ println("1. Autonomous System (No Forcing)")
println("-"^60)
# Standard Van der Pol with moderate damping
-model = VanDerPol(ϵ=1.0, F=0.0, ω=0.0)
+model = VanDerPol(; ϵ = 1.0, F = 0.0, ω = 0.0)
println("\nModel parameters:")
println(" ϵ (damping) = $(model.ϵ)")
println(" F (forcing amplitude) = $(model.F)")
@@ -29,7 +29,7 @@ println("\nInitial condition: ", x0)
# Lyapunov exponent for limit cycle
println("\nCalculating Lyapunov exponent...")
-λ_max = lyapunov_exponent(model, x0, 0.1; warmup=1000, n_iterations=3000)
+λ_max = lyapunov_exponent(model, x0, 0.1; warmup = 1000, n_iterations = 3000)
println("Largest Lyapunov exponent: ", λ_max)
println("\nInterpretation:")
@@ -52,9 +52,9 @@ println()
ϵ_values = [0.1, 1.0, 5.0, 10.0]
for ϵ in ϵ_values
- test_model = VanDerPol(ϵ=ϵ)
- λ = lyapunov_exponent(test_model, x0, 0.1; warmup=500, n_iterations=1000)
-
+ test_model = VanDerPol(; ϵ = ϵ)
+ λ = lyapunov_exponent(test_model, x0, 0.1; warmup = 500, n_iterations = 1000)
+
oscillation_type = if ϵ < 0.5
"Nearly sinusoidal"
elseif ϵ < 3.0
@@ -62,8 +62,8 @@ for ϵ in ϵ_values
else
"Relaxation oscillations"
end
-
- println(" ϵ = ", ϵ, "\tλ ≈ ", round(λ, digits=4), "\t[", oscillation_type, "]")
+
+ println(" ϵ = ", ϵ, "\tλ ≈ ", round(λ; digits = 4), "\t[", oscillation_type, "]")
end
println("\nObservations:")
@@ -76,7 +76,7 @@ println("3. Forced Van der Pol Oscillator")
println("-"^60)
# Forced oscillator
-forced_model = VanDerPol(ϵ=1.0, F=5.0, ω=2.466)
+forced_model = VanDerPol(; ϵ = 1.0, F = 5.0, ω = 2.466)
println("\nForced oscillator parameters:")
println(" ϵ = $(forced_model.ϵ)")
println(" F = $(forced_model.F)")
diff --git a/script/makefigures.jl b/script/makefigures.jl
index 4e936ad..13ce135 100644
--- a/script/makefigures.jl
+++ b/script/makefigures.jl
@@ -4,7 +4,11 @@ using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using DynamicalModels
-using Plots, LinearAlgebra, Random, Statistics, FFTW
+using FFTW
+using LinearAlgebra
+using Plots
+using Random
+using Statistics
const FIG_DIR = joinpath(pkgdir(DynamicalModels), "docs", "src", "assets", "figures")
const FIG_ODE = joinpath(FIG_DIR, "solver", "ode_solver")
@@ -29,4 +33,4 @@ println("Passed arguments ARGS = $(test_args) to tests.")
include(filepath)
end
end
-end
\ No newline at end of file
+end
diff --git a/script/model/fig_henonmap.jl b/script/model/fig_henonmap.jl
index 3046e3a..5512f47 100644
--- a/script/model/fig_henonmap.jl
+++ b/script/model/fig_henonmap.jl
@@ -10,25 +10,71 @@ let
ncols = 2
nrows = 2
- plt = plot(layout=(2, 2), size=(300 * ncols, 300 * nrows), xlabel="", ylabel="", legend=false)
- scatter!(plt, x_sol[:, 1], x_sol[:, 2], xlim=(-1.5, 1.5), ylim=(-0.5, 0.5), markersize=0.5, subplot=1)
- scatter!(plt, x_sol[:, 1], x_sol[:, 2], xlim=(0.0, 0.6), ylim=(0.1, 0.3), markersize=0.5, subplot=2)
- scatter!(plt, x_sol[:, 1], x_sol[:, 2], xlim=(0.1, 0.2), ylim=(0.2, 0.25), markersize=0.5, subplot=3)
- scatter!(plt, x_sol[:, 1], x_sol[:, 2], xlim=(0.105, 0.12), ylim=(0.23, 0.245), markersize=0.5, subplot=4)
+ plt = plot(;
+ layout = (2, 2),
+ size = (300 * ncols, 300 * nrows),
+ xlabel = "",
+ ylabel = "",
+ legend = false
+ )
+ scatter!(
+ plt,
+ x_sol[:, 1],
+ x_sol[:, 2];
+ xlim = (-1.5, 1.5),
+ ylim = (-0.5, 0.5),
+ markersize = 0.5,
+ subplot = 1
+ )
+ scatter!(
+ plt,
+ x_sol[:, 1],
+ x_sol[:, 2];
+ xlim = (0.0, 0.6),
+ ylim = (0.1, 0.3),
+ markersize = 0.5,
+ subplot = 2
+ )
+ scatter!(
+ plt,
+ x_sol[:, 1],
+ x_sol[:, 2];
+ xlim = (0.1, 0.2),
+ ylim = (0.2, 0.25),
+ markersize = 0.5,
+ subplot = 3
+ )
+ scatter!(
+ plt,
+ x_sol[:, 1],
+ x_sol[:, 2];
+ xlim = (0.105, 0.12),
+ ylim = (0.23, 0.245),
+ markersize = 0.5,
+ subplot = 4
+ )
savefig(plt, joinpath(FIG_HENON, "00_henon_map.png"))
# TODO: Divergence of nearby trajectories
x0 = randn(2)
x1 = randn(2)
- diff_list = [1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3]
- plt = plot(xlabel="Steps", ylabel="Distance", title="Divergence of nearby trajectories in Henon Map", yscale=:log10, label="", legend=:topleft, minorgrid=true)
+ diff_list = [1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3]
+ plt = plot(;
+ xlabel = "Steps",
+ ylabel = "Distance",
+ title = "Divergence of nearby trajectories in Henon Map",
+ yscale = :log10,
+ label = "",
+ legend = :topleft,
+ minorgrid = true
+ )
for dd in diff_list
xx = x0 .+ x1 .* dd
n_steps = 50
x_sol0 = map_solver(model, n_steps, x0)
x_sol1 = map_solver(model, n_steps, xx)
norms = [norm(x_sol1[n, :] .- x_sol0[n, :]) for n in 1:size(x_sol0, 1)]
- plot!(plt, 1:n_steps, norms; label="initial diff $(dd)", lw=0.5)
+ plot!(plt, 1:n_steps, norms; label = "initial diff $(dd)", lw = 0.5)
end
savefig(plt, joinpath(FIG_HENON, "01_henon_divergence_initial_diff.png"))
-end
\ No newline at end of file
+end
diff --git a/script/model/fig_lorenz.jl b/script/model/fig_lorenz.jl
index d6d7abc..76fce94 100644
--- a/script/model/fig_lorenz.jl
+++ b/script/model/fig_lorenz.jl
@@ -10,23 +10,34 @@ let
x_sol = DynamicalModels.ode_solver(DynamicalModels.RK4, model, t_list, x0)
filter = t_list .> 5
x_sol_plot = @view x_sol[filter, :]
- plt = plot3d(x_sol_plot[:, 1], x_sol_plot[:, 2], x_sol_plot[:, 3];
- array=true, xlabel="X", ylabel="Y", zlabel="Z",
- title="Attractor of " * string(typeof(model)), aspect_ratio=:equal, label="", lw=0.50)
+ plt = plot3d(
+ x_sol_plot[:, 1], x_sol_plot[:, 2], x_sol_plot[:, 3];
+ array = true, xlabel = "X", ylabel = "Y", zlabel = "Z",
+ title = "Attractor of " * string(typeof(model)), aspect_ratio = :equal,
+ label = "", lw = 0.5
+ )
savefig(plt, joinpath(FIG_LORENZ, "00_lorenz_attractor.png"))
# TODO : Divergence of nearby trajectories
- diff_list = [1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3]
+ diff_list = [1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3]
x1 = [randn(), randn(), randn()]
x_sol = DynamicalModels.ode_solver(DynamicalModels.RK4, model, t_list, x0)
- plt = plot(xlabel="Time", ylabel="Distance", title="Divergence of nearby trajectories in " * string(typeof(model)), yscale=:log10, label="", minorgrid=true, legend=:topleft)
+ plt = plot(;
+ xlabel = "Time",
+ ylabel = "Distance",
+ title = "Divergence of nearby trajectories in " * string(typeof(model)),
+ yscale = :log10,
+ label = "",
+ minorgrid = true,
+ legend = :topleft
+ )
for dd in diff_list
xx = x0 .+ x1 .* dd
x_sol1 = DynamicalModels.ode_solver(DynamicalModels.RK4, model, t_list, xx)
diff = x_sol1 .- x_sol
norms = [norm(diff[i, :]) for i in 1:size(diff, 1)]
- plot!(t_list, norms; label="initial diff $(dd)", lw=0.5)
+ plot!(t_list, norms; label = "initial diff $(dd)", lw = 0.5)
end
- savefig(plt, joinpath(FIG_LORENZ, "01_lorenz_divergence_initial_diff.png"))
-end
\ No newline at end of file
+ savefig(plt, joinpath(FIG_LORENZ, "01_lorenz_divergence_initial_diff.png"))
+end
diff --git a/script/model/fig_rossler.jl b/script/model/fig_rossler.jl
index 0f44c8c..5e25a6c 100644
--- a/script/model/fig_rossler.jl
+++ b/script/model/fig_rossler.jl
@@ -10,20 +10,40 @@ let
x_sol = DynamicalModels.ode_solver(DynamicalModels.RK4, model, t_list, x0)
filter = t_list .> 5
x_sol_plot = @view x_sol[filter, :]
- plt = plot3d(x_sol_plot[:, 1], x_sol_plot[:, 2], x_sol_plot[:, 3]; array=true, xlabel="X", ylabel="Y", zlabel="Z", title="Attractor of " * string(typeof(model)), aspect_ratio=:equal, label="", lw=0.50)
+ plt = plot3d(
+ x_sol_plot[:, 1],
+ x_sol_plot[:, 2],
+ x_sol_plot[:, 3];
+ array = true,
+ xlabel = "X",
+ ylabel = "Y",
+ zlabel = "Z",
+ title = "Attractor of " * string(typeof(model)),
+ aspect_ratio = :equal,
+ label = "",
+ lw = 0.5
+ )
savefig(plt, joinpath(FIG_ROSSLER, "00_rossler_attractor.png"))
# TODO : Divergence of nearby trajectories
- diff_list = [1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3]
+ diff_list = [1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3]
x1 = [randn(), randn(), randn()]
x_sol = DynamicalModels.ode_solver(DynamicalModels.RK4, model, t_list, x0)
- plt = plot(xlabel="Time", ylabel="Distance", title="Divergence of nearby trajectories in " * string(typeof(model)), yscale=:log10, label="", legend=:topleft, minorgrid=true)
+ plt = plot(;
+ xlabel = "Time",
+ ylabel = "Distance",
+ title = "Divergence of nearby trajectories in " * string(typeof(model)),
+ yscale = :log10,
+ label = "",
+ legend = :topleft,
+ minorgrid = true
+ )
for dd in diff_list
xx = x0 .+ x1 .* dd
x_sol1 = DynamicalModels.ode_solver(DynamicalModels.RK4, model, t_list, xx)
diff = x_sol1 .- x_sol
norms = [norm(diff[i, :]) for i in 1:size(diff, 1)]
- plot!(t_list, norms; label="initial diff $(dd)", lw=0.5)
+ plot!(t_list, norms; label = "initial diff $(dd)", lw = 0.5)
end
savefig(plt, joinpath(FIG_ROSSLER, "01_rossler_divergence_initial_diff.png"))
-end
\ No newline at end of file
+end
diff --git a/script/model/fig_vanderpol.jl b/script/model/fig_vanderpol.jl
index 3a86140..ce370ab 100644
--- a/script/model/fig_vanderpol.jl
+++ b/script/model/fig_vanderpol.jl
@@ -1,6 +1,6 @@
const FIG_VANDERPOL = joinpath(FIG_MODEL, "vanderpol")
mkpath(FIG_VANDERPOL)
-let
+let
model = VanDerPol()
x0 = [randn(), randn()]
dt = 0.01
@@ -12,19 +12,19 @@ let
filter = t_list .> 5
ϵ_list = [0.0, 0.1, 1.0, 2.0, 4.0]
- plt1 = plot(title="Van der Pol Oscillator", xlabel="Time", ylabel="x",)
- plt2 = plot(title="Van der Pol Oscillator Phase Space", xlabel="x", ylabel="v",)
- plt3 = plot(title="Van der Pol Poincare Section", xlabel="x_n", ylabel="x_n+1",)
- plt4 = plot(title="Van der Pol return maps", xlabel="x_n", ylabel="x_n+1",)
- plt5 = plot(title="Amplitudes", xlabel="ϵ", ylabel="Amplitude",)
- plt6 = plot(title="Amplitudes", xlabel="ϵ", ylabel="Amplitude",)
+ plt1 = plot(; title = "Van der Pol Oscillator", xlabel = "Time", ylabel = "x")
+ plt2 = plot(; title = "Van der Pol Oscillator Phase Space", xlabel = "x", ylabel = "v")
+ plt3 = plot(; title = "Van der Pol Poincare Section", xlabel = "x_n", ylabel = "x_n+1")
+ plt4 = plot(; title = "Van der Pol return maps", xlabel = "x_n", ylabel = "x_n+1")
+ plt5 = plot(; title = "Amplitudes", xlabel = "ϵ", ylabel = "Amplitude")
+ plt6 = plot(; title = "Amplitudes", xlabel = "ϵ", ylabel = "Amplitude")
function poincare_section(x_sol, t_list)
Poincare = []
for (t, tt) in enumerate(t_list)
if t == length(t_list)
break
end
- if x_sol[t, 2] * x_sol[t+1, 2] < 0 && x_sol[t, 1] > 0
+ if x_sol[t, 2] * x_sol[t + 1, 2] < 0 && x_sol[t, 1] > 0
push!(Poincare, x_sol[t, 1])
end
end
@@ -32,7 +32,7 @@ let
end
for ϵ in ϵ_list
# TODO solve
- model = VanDerPol(ϵ=ϵ)
+ model = VanDerPol(; ϵ = ϵ)
x0 = [2.0, 0.0]
t_max = 20.0
dt = 0.01
@@ -41,8 +41,8 @@ let
# TODO : solve nulclines
# plot!(plt2, ~~~)
# TODO : show plots
- plot!(plt1, t_list, x_sol[:, 1], label="ϵ=$ϵ")
- plot!(plt2, x_sol[:, 1], x_sol[:, 2], label="ϵ=$ϵ")
+ plot!(plt1, t_list, x_sol[:, 1]; label = "ϵ=$ϵ")
+ plot!(plt2, x_sol[:, 1], x_sol[:, 2]; label = "ϵ=$ϵ")
# TODO : Poincare 断面
x0 = [randn(), randn()]
t_max = 200.0
@@ -50,24 +50,24 @@ let
t_list = collect(0:dt:t_max)
x_sol = DynamicalModels.ode_solver(DynamicalModels.RK4, model, t_list, x0)
Poincare = poincare_section(x_sol, t_list)
- plot!(plt3, Poincare[1:end-1], Poincare[2:end], label="ϵ=$ϵ")
+ plot!(plt3, Poincare[1:(end - 1)], Poincare[2:end]; label = "ϵ=$ϵ")
end
ϵ_list = collect(0.0:0.1:4.0)
for ϵ in ϵ_list
- model = VanDerPol(ϵ=ϵ)
+ model = VanDerPol(; ϵ = ϵ)
x0 = [randn(), randn()]
t_max = 200.0
dt = 0.01
t_list = collect(0:dt:t_max)
x_sol = DynamicalModels.ode_solver(DynamicalModels.RK4, model, t_list, x0)
Poincare = poincare_section(x_sol, t_list)
- plot!(plt4, Poincare[1:end-1], Poincare[2:end], label="")
+ plot!(plt4, Poincare[1:(end - 1)], Poincare[2:end]; label = "")
end
- plot!(plt4, [0, 3], [0, 3], label="y=x", linestyle=:dash, color=:black)
+ plot!(plt4, [0, 3], [0, 3]; label = "y=x", linestyle = :dash, color = :black)
ϵ_list = collect(0.01:0.01:2.0)
Amplitude = zeros(length(ϵ_list))
for (i, ϵ) in enumerate(ϵ_list)
- model = VanDerPol(ϵ=ϵ)
+ model = VanDerPol(; ϵ = ϵ)
x0 = [randn(), randn()]
t_max = 200.0
dt = 0.01
@@ -76,9 +76,16 @@ let
filter = t_list .> (t_max / 2)
Amplitude[i] = maximum(x_sol[filter, 1]) - minimum(x_sol[filter, 1])
end
- plot!(plt5, ϵ_list, Amplitude, label="Amplitude")
- plot!(plt6, ϵ_list, Amplitude, label="Amplitude", xscale=:log10, yscale=:log10)
- plot!(plt6, ϵ_list, sqrt.(ϵ_list) * Amplitude[end] / sqrt(ϵ_list[end]), label="ϵ^0.5", xscale=:log10, yscale=:log10)
+ plot!(plt5, ϵ_list, Amplitude; label = "Amplitude")
+ plot!(plt6, ϵ_list, Amplitude; label = "Amplitude", xscale = :log10, yscale = :log10)
+ plot!(
+ plt6,
+ ϵ_list,
+ sqrt.(ϵ_list) * Amplitude[end] / sqrt(ϵ_list[end]);
+ label = "ϵ^0.5",
+ xscale = :log10,
+ yscale = :log10
+ )
savefig(plt1, joinpath(FIG_VANDERPOL, "00_vanderpol_time.png"))
savefig(plt2, joinpath(FIG_VANDERPOL, "01_vanderpol_phase_space.png"))
savefig(plt3, joinpath(FIG_VANDERPOL, "02_vanderpol_poincare_section.png"))
@@ -92,14 +99,32 @@ let
for ϵ in ϵ_list
ncols = length(F_list)
nrows = length(ω_list)
- plt1 = plot(layout=(nrows, ncols), size=(300 * ncols, 250 * nrows), xlabel="", ylabel="", legend=false)
- plt2 = plot(layout=(nrows, ncols), size=(300 * ncols, 250 * nrows), xlabel="", ylabel="", legend=false)
- plt3 = plot(layout=(nrows, ncols), size=(300 * ncols, 250 * nrows), xlabel="", ylabel="", legend=false)
+ plt1 = plot(;
+ layout = (nrows, ncols),
+ size = (300 * ncols, 250 * nrows),
+ xlabel = "",
+ ylabel = "",
+ legend = false
+ )
+ plt2 = plot(;
+ layout = (nrows, ncols),
+ size = (300 * ncols, 250 * nrows),
+ xlabel = "",
+ ylabel = "",
+ legend = false
+ )
+ plt3 = plot(;
+ layout = (nrows, ncols),
+ size = (300 * ncols, 250 * nrows),
+ xlabel = "",
+ ylabel = "",
+ legend = false
+ )
for (f, F) in enumerate(F_list)
for (w, ω) in enumerate(ω_list)
idx = (f - 1) * ncols + w
# solve
- model = VanDerPol(ϵ=ϵ, F=F, ω=ω)
+ model = VanDerPol(; ϵ = ϵ, F = F, ω = ω)
x0 = [1.0, 0.0]
t_max = 200.0
dt = 0.04
@@ -109,16 +134,30 @@ let
filter = t_list .> (t_max / 2)
xs = x_sol[filter, 1]
ys = x_sol[filter, 2]
- plot!(plt1, xs, ys; subplot=idx, lw=1.5)
- plot!(plt1, title="F=$(F), ω=$(ω)", subplot=idx)
+ plot!(plt1, xs, ys; subplot = idx, lw = 1.5)
+ plot!(plt1; title = "F=$(F), ω=$(ω)", subplot = idx)
# time position
filter = t_list .> (t_max - 40.0)
xs = x_sol[filter, 1]
ys = x_sol[filter, 2]
- plot!(plt2, t_list[filter], xs; title="F=$(F), ω=$(ω)", subplot=idx, lw=1.5)
- plot!(plt2, t_list[filter], ys; title="F=$(F), ω=$(ω)", subplot=idx, lw=1.5)
+ plot!(
+ plt2,
+ t_list[filter],
+ xs;
+ title = "F=$(F), ω=$(ω)",
+ subplot = idx,
+ lw = 1.5
+ )
+ plot!(
+ plt2,
+ t_list[filter],
+ ys;
+ title = "F=$(F), ω=$(ω)",
+ subplot = idx,
+ lw = 1.5
+ )
# FFT
- model = VanDerPol(ϵ=ϵ, F=F, ω=ω)
+ model = VanDerPol(; ϵ = ϵ, F = F, ω = ω)
x0 = [1.0, 0.0]
t_max = 4000.0
dt = 0.04
@@ -129,8 +168,8 @@ let
xf = real.(fft(x_sol[:, 1]))
freq = fftshift(freq)
xf = fftshift(xf)
- plot!(plt3, freq[freq.>=0], abs.(xf[freq.>=0]); subplot=idx, lw=1.5)
- plot!(plt3, title="F=$(F), ω=$(ω)", subplot=idx, xlim=(0, 4.0))
+ plot!(plt3, freq[freq .>= 0], abs.(xf[freq .>= 0]); subplot = idx, lw = 1.5)
+ plot!(plt3; title = "F=$(F), ω=$(ω)", subplot = idx, xlim = (0, 4.0))
end
end
savefig(plt1, joinpath(FIG_VANDERPOL, "10_vanderpol_phase_space_ϵ$(ϵ).png"))
@@ -149,7 +188,7 @@ let
arnoldi = zeros(nF, nω)
for (f, F) in enumerate(F_list)
for (w, ω) in enumerate(ω_list)
- model = VanDerPol(ϵ=ϵ, F=F, ω=ω)
+ model = VanDerPol(; ϵ = ϵ, F = F, ω = ω)
x0 = [1.0, 0.0]
t_max = 1000.0
dt = 0.1
@@ -164,16 +203,25 @@ let
if i == length(tt)
break
end
- if xs[i, 2] * xs[i+1, 2] < 0 && xs[i, 1] > 0
+ if xs[i, 2] * xs[i + 1, 2] < 0 && xs[i, 1] > 0
push!(times, t)
end
end
- time_diff = times[2:end] .- times[1:end-1]
+ time_diff = times[2:end] .- times[1:(end - 1)]
arnoldi[f, w] = 2π / mean(time_diff)
arnoldi[f, w] /= ω
end
end
- plt = heatmap(ω_list, F_list, arnoldi; ylabel="F", xlabel="ω", title="Arnold Tongue ϵ=$(ϵ)", colorbar_title="relative frequency", clim=(0, 3))
+ plt = heatmap(
+ ω_list,
+ F_list,
+ arnoldi;
+ ylabel = "F",
+ xlabel = "ω",
+ title = "Arnold Tongue ϵ=$(ϵ)",
+ colorbar_title = "relative frequency",
+ clim = (0, 3)
+ )
savefig(plt, joinpath(FIG_VANDERPOL, "20_vanderpol_arnold_tongue_ϵ$(ϵ).png"))
end
-end
\ No newline at end of file
+end
diff --git a/script/solver/fig_odesolver.jl b/script/solver/fig_odesolver.jl
index 6adb3e6..564cd8e 100644
--- a/script/solver/fig_odesolver.jl
+++ b/script/solver/fig_odesolver.jl
@@ -1,5 +1,5 @@
let
- model = HarmonicOscillator(k=1.0, m=1.0)
+ model = HarmonicOscillator(; k = 1.0, m = 1.0)
x0 = [1.0, 0.0]
t_max = 10.0
dt = 0.01
@@ -7,24 +7,38 @@ let
# TODO : calculate convergence order and plot
N_list = [20, 40, 80, 160, 320, 640, 1280]
- plt = plot(
- title="ODE Solver Convergence Order",
- xaxis=:log, yaxis=:log,
- xlabel="Step Size (h)",
- ylabel="Final Error (E)",
- legend=:bottomright,
- minorgrid=true
+ plt = plot(;
+ title = "ODE Solver Convergence Order",
+ xaxis = :log, yaxis = :log,
+ xlabel = "Step Size (h)",
+ ylabel = "Final Error (E)",
+ legend = :bottomright,
+ minorgrid = true
+ )
+ h_ref = [1.0e-3, 1.0]
+ plot!(plt, h_ref, h_ref .^ 1; linestyle = :dash, color = :gray, label = "Order 1 (h)")
+ plot!(
+ plt,
+ h_ref,
+ h_ref .^ 2;
+ linestyle = :dash,
+ color = :darkgray,
+ label = "Order 2 (h^2)"
+ )
+ plot!(
+ plt,
+ h_ref,
+ h_ref .^ 4;
+ linestyle = :dash,
+ color = :black,
+ label = "Order 4 (h^4)"
)
- h_ref = [1e-3, 1.0]
- plot!(plt, h_ref, h_ref .^ 1, linestyle=:dash, color=:gray, label="Order 1 (h)")
- plot!(plt, h_ref, h_ref .^ 2, linestyle=:dash, color=:darkgray, label="Order 2 (h^2)")
- plot!(plt, h_ref, h_ref .^ 4, linestyle=:dash, color=:black, label="Order 4 (h^4)")
for solver_func in DynamicalModels.ODE_SOLVERS
solver_name = string(solver_func)
h_list = Float64[]
error_list = Float64[]
for N in N_list
- t_list = range(0, t_max; length=N)
+ t_list = range(0, t_max; length = N)
h = t_list[2] - t_list[1]
x_list_num = DynamicalModels.ode_solver(solver_func, model, t_list, x0)
sol_exact_list = solve_exact(model, t_list, x0)
@@ -33,27 +47,34 @@ let
push!(h_list, h)
push!(error_list, final_error)
end
- plot!(plt, h_list, error_list, marker=:o, label=solver_name)
+ plot!(plt, h_list, error_list; marker = :o, label = solver_name)
end
savefig(plt, joinpath(FIG_ODE, "convergence_order.png"))
# TODO : plot trajectories for different solvers
- plt_traj = plot(
- title="ODE Solver Trajectories",
- xlabel="Time (t)",
- ylabel="Position (x)",
- legend=:topright,
- minorgrid=true
+ plt_traj = plot(;
+ title = "ODE Solver Trajectories",
+ xlabel = "Time (t)",
+ ylabel = "Position (x)",
+ legend = :topright,
+ minorgrid = true
)
N = 1000
- t_list = range(0, t_max; length=N)
+ t_list = range(0, t_max; length = N)
for solver_func in DynamicalModels.ODE_SOLVERS
solver_name = string(solver_func)
x_list_num = DynamicalModels.ode_solver(solver_func, model, t_list, x0)
- plot!(plt_traj, t_list, x_list_num[:, 1], label=solver_name)
+ plot!(plt_traj, t_list, x_list_num[:, 1]; label = solver_name)
end
sol_exact_list = solve_exact(model, t_list, x0)
x_list_exact = hcat(sol_exact_list...)' # [N, 2]
- plot!(plt_traj, t_list, x_list_exact[:, 1], linestyle=:dash, color=:black, label="Exact Solution")
+ plot!(
+ plt_traj,
+ t_list,
+ x_list_exact[:, 1];
+ linestyle = :dash,
+ color = :black,
+ label = "Exact Solution"
+ )
savefig(plt_traj, joinpath(FIG_ODE, "trajectories_harmonic_oscillator.png"))
-end
\ No newline at end of file
+end
diff --git a/src/abstractdynamicalsystem.jl b/src/abstractdynamicalsystem.jl
index 40657c3..4788fd6 100644
--- a/src/abstractdynamicalsystem.jl
+++ b/src/abstractdynamicalsystem.jl
@@ -31,6 +31,7 @@ export AbstractMap
"""
AbstractStochasticSystem{T} <: AbstractSystem{T}
+
Abstract type for stochastic systems. It is intended to be used for models
that employ stochastic methods such as the Metropolis algorithm or Gibbs sampling.
"""
@@ -39,6 +40,7 @@ export AbstractStochasticSystem
"""
AbstractSDEModel{T} <: AbstractStochasticSystem{T}
+
Abstract type for stochastic differential equation models. It is intended to be passed to `sde_solver()`.
For `params :: AbstractSDEModel{T}`, a stochastic differential equation `params(t, x) -> (drift, diffusion)` is defined,
where `(t, x)` are the variables.
diff --git a/src/analysis/dimension.jl b/src/analysis/dimension.jl
index ecfd65e..7aee007 100644
--- a/src/analysis/dimension.jl
+++ b/src/analysis/dimension.jl
@@ -1,4 +1,4 @@
-const NUMERICAL_EPSILON = 1e-10
+const NUMERICAL_EPSILON = 1.0e-10
"""
correlation_dimension(points; r_min=0.01, r_max=10.0, n_r=50)
@@ -6,41 +6,46 @@ const NUMERICAL_EPSILON = 1e-10
Estimate the correlation dimension of a set of points using the Grassberger-Procaccia algorithm.
# Arguments
-- `points`: Vector of points (each point is a vector)
-- `r_min`: Minimum radius for correlation sum (default: 0.01)
-- `r_max`: Maximum radius for correlation sum (default: 10.0)
-- `n_r`: Number of radius values to test (default: 50)
+
+ - `points`: Vector of points (each point is a vector)
+ - `r_min`: Minimum radius for correlation sum (default: 0.01)
+ - `r_max`: Maximum radius for correlation sum (default: 10.0)
+ - `n_r`: Number of radius values to test (default: 50)
# Returns
-- Tuple of (radii, correlation_sums, estimated_dimension)
+
+ - Tuple of (radii, correlation_sums, estimated_dimension)
The correlation dimension is estimated from the slope of log(C(r)) vs log(r).
# Example
+
```julia
model = Lorenz()
trajectory = solve_trajectory(model, [1.0, 1.0, 1.0], 1000.0, 0.01)
radii, C, dim = correlation_dimension(trajectory)
```
"""
-function correlation_dimension(points::Vector{Vector{T}};
- r_min=0.01, r_max=10.0, n_r=50) where T
+function correlation_dimension(
+ points::Vector{Vector{T}};
+ r_min = 0.01, r_max = 10.0, n_r = 50
+ ) where {T}
N = length(points)
-
+
if N < 2
error("Need at least 2 points to calculate correlation dimension")
end
-
+
# Generate logarithmically spaced radii
- radii = exp.(range(log(r_min), log(r_max), length=n_r))
-
+ radii = exp.(range(log(r_min), log(r_max); length = n_r))
+
# Calculate correlation sums
C = zeros(Float64, n_r)
-
+
for (i, r) in enumerate(radii)
count = 0
for j in 1:N
- for k in (j+1):N
+ for k in (j + 1):N
if norm(points[j] - points[k]) < r
count += 1
end
@@ -49,38 +54,38 @@ function correlation_dimension(points::Vector{Vector{T}};
# Correlation sum
C[i] = 2.0 * count / (N * (N - 1))
end
-
+
# Estimate dimension from linear region
# Use middle portion of the curve to avoid edge effects
log_r = log.(radii)
log_C = log.(C .+ NUMERICAL_EPSILON) # Add small value to avoid log(0)
-
+
# Find linear region (middle 60% of data where C > 0)
valid_indices = findall(C .> NUMERICAL_EPSILON)
if length(valid_indices) < 2
return radii, C, NaN
end
-
+
mid_start = valid_indices[div(length(valid_indices), 5)]
mid_end = valid_indices[div(4 * length(valid_indices), 5)]
-
+
if mid_end <= mid_start
return radii, C, NaN
end
-
+
# Linear regression on log-log plot
X = log_r[mid_start:mid_end]
Y = log_C[mid_start:mid_end]
-
+
# Calculate slope (dimension estimate)
n_points = length(X)
sum_x = sum(X)
sum_y = sum(Y)
sum_xy = sum(X .* Y)
sum_x2 = sum(X .^ 2)
-
+
dimension = (n_points * sum_xy - sum_x * sum_y) / (n_points * sum_x2 - sum_x^2)
-
+
return radii, C, dimension
end
export correlation_dimension
@@ -91,10 +96,12 @@ export correlation_dimension
Calculate the Kaplan-Yorke (Lyapunov) dimension from a spectrum of Lyapunov exponents.
# Arguments
-- `lyapunov_exponents`: Vector of Lyapunov exponents (should be sorted in descending order)
+
+ - `lyapunov_exponents`: Vector of Lyapunov exponents (should be sorted in descending order)
# Returns
-- Kaplan-Yorke dimension
+
+ - Kaplan-Yorke dimension
The Kaplan-Yorke dimension is defined as:
D_KY = j + (λ_1 + λ_2 + ... + λ_j) / |λ_(j+1)|
@@ -102,34 +109,35 @@ D_KY = j + (λ_1 + λ_2 + ... + λ_j) / |λ_(j+1)|
where j is the largest integer such that λ_1 + λ_2 + ... + λ_j ≥ 0.
# Example
+
```julia
model = Lorenz()
x0 = [1.0, 1.0, 1.0]
λs = lyapunov_spectrum(model, x0, 0.1)
D_KY = kaplan_yorke_dimension(λs)
-```
+``` # Ensure sorted in descending order
"""
-function kaplan_yorke_dimension(lyapunov_exponents::Vector{T}) where T
+function kaplan_yorke_dimension(lyapunov_exponents::Vector{T}) where {T}
# Ensure sorted in descending order
- λs = sort(lyapunov_exponents, rev=true)
-
+ λs = sort(lyapunov_exponents; rev = true)
+
# Find largest j such that sum of first j exponents is non-negative
cumsum_λ = cumsum(λs)
j = findlast(cumsum_λ .>= 0)
-
+
if isnothing(j)
# All exponents are negative
return 0.0
end
-
+
if j == length(λs)
# All exponents sum to positive (unbounded growth)
return Float64(length(λs))
end
-
+
# Calculate Kaplan-Yorke dimension
- D_KY = j + cumsum_λ[j] / abs(λs[j+1])
-
+ D_KY = j + cumsum_λ[j] / abs(λs[j + 1])
+
return D_KY
end
export kaplan_yorke_dimension
@@ -140,65 +148,70 @@ export kaplan_yorke_dimension
Estimate the box-counting (capacity) dimension of a trajectory.
# Arguments
-- `trajectory`: Vector of points representing a trajectory
-- `box_sizes`: Vector of box sizes to test (if not provided, automatically generated)
+
+ - `trajectory`: Vector of points representing a trajectory
+ - `box_sizes`: Vector of box sizes to test (if not provided, automatically generated)
# Returns
-- Tuple of (box_sizes, N_boxes, estimated_dimension)
+
+ - Tuple of (box_sizes, N_boxes, estimated_dimension)
# Example
+
```julia
model = Lorenz()
trajectory = solve_trajectory(model, [1.0, 1.0, 1.0], 1000.0, 0.01)
sizes, N, dim = box_counting_dimension(trajectory)
```
"""
-function box_counting_dimension(trajectory::Vector{Vector{T}},
- box_sizes=nothing) where T
+function box_counting_dimension(
+ trajectory::Vector{Vector{T}},
+ box_sizes = nothing
+ ) where {T}
if isempty(trajectory)
error("Trajectory is empty")
end
-
+
dim = length(trajectory[1])
-
+
# Find bounding box
- mins = minimum(hcat(trajectory...), dims=2)[:]
- maxs = maximum(hcat(trajectory...), dims=2)[:]
+ mins = minimum(hcat(trajectory...); dims = 2)[:]
+ maxs = maximum(hcat(trajectory...); dims = 2)[:]
extent = maximum(maxs - mins)
-
+
# Generate box sizes if not provided
if isnothing(box_sizes)
box_sizes = extent ./ (2 .^ (1:10))
end
-
+
N_boxes = zeros(Int, length(box_sizes))
-
+
for (i, ε) in enumerate(box_sizes)
# Count number of occupied boxes
occupied = Set{Vector{Int}}()
-
+
for point in trajectory
# Determine which box this point falls into
box_indices = floor.(Int, (point - mins) / ε)
push!(occupied, box_indices)
end
-
+
N_boxes[i] = length(occupied)
end
-
+
# Estimate dimension from slope of log(N) vs log(1/ε)
log_inv_ε = log.(1 ./ box_sizes)
log_N = log.(Float64.(N_boxes))
-
+
# Linear regression
n_points = length(log_inv_ε)
sum_x = sum(log_inv_ε)
sum_y = sum(log_N)
sum_xy = sum(log_inv_ε .* log_N)
sum_x2 = sum(log_inv_ε .^ 2)
-
+
dimension = (n_points * sum_xy - sum_x * sum_y) / (n_points * sum_x2 - sum_x^2)
-
+
return box_sizes, N_boxes, dimension
end
export box_counting_dimension
diff --git a/src/analysis/lyapunov.jl b/src/analysis/lyapunov.jl
index d4a558e..33c6321 100644
--- a/src/analysis/lyapunov.jl
+++ b/src/analysis/lyapunov.jl
@@ -1,4 +1,4 @@
-const PERTURBATION_SIZE = 1e-8
+const PERTURBATION_SIZE = 1.0e-8
"""
lyapunov_exponent(model, x0, time_step; dt=0.01, warmup=1000, n_iterations=10000)
@@ -6,30 +6,35 @@ const PERTURBATION_SIZE = 1e-8
Calculate the largest Lyapunov exponent for a continuous dynamical system.
# Arguments
-- `model`: The dynamical model (must be callable with signature `model(t, x)`)
-- `x0`: Initial condition vector
-- `time_step`: Total integration time per iteration
-- `dt`: Time step for integration (default: 0.01)
-- `warmup`: Number of warmup iterations before calculating exponent (default: 1000)
-- `n_iterations`: Number of iterations for calculation (default: 10000)
+
+ - `model`: The dynamical model (must be callable with signature `model(t, x)`)
+ - `x0`: Initial condition vector
+ - `time_step`: Total integration time per iteration
+ - `dt`: Time step for integration (default: 0.01)
+ - `warmup`: Number of warmup iterations before calculating exponent (default: 1000)
+ - `n_iterations`: Number of iterations for calculation (default: 10000)
# Returns
-- Largest Lyapunov exponent
+
+ - Largest Lyapunov exponent
# Example
+
```julia
model = Lorenz()
x0 = [1.0, 1.0, 1.0]
λ = lyapunov_exponent(model, x0, 0.1)
```
"""
-function lyapunov_exponent(model, x0::AbstractVector{T}, time_step::Real;
- dt=0.01, warmup=1000, n_iterations=10000) where T
+function lyapunov_exponent(
+ model, x0::AbstractVector{T}, time_step::Real;
+ dt = 0.01, warmup = 1000, n_iterations = 10000
+ ) where {T}
dim = length(x0)
x = copy(x0)
δx = randn(T, dim)
δx = δx / norm(δx) * PERTURBATION_SIZE # Small perturbation
-
+
# Warmup phase
for _ in 1:warmup
x, _ = rk4_step(model, 0.0, x, time_step, dt)
@@ -37,21 +42,21 @@ function lyapunov_exponent(model, x0::AbstractVector{T}, time_step::Real;
δx = δx_evolved - x
δx = δx / norm(δx) * PERTURBATION_SIZE
end
-
+
# Calculation phase
sum_log = 0.0
for _ in 1:n_iterations
x_new, _ = rk4_step(model, 0.0, x, time_step, dt)
δx_evolved, _ = rk4_step(model, 0.0, x + δx, time_step, dt)
δx_new = δx_evolved - x_new
-
+
d = norm(δx_new)
sum_log += log(d / PERTURBATION_SIZE)
-
+
x = x_new
δx = δx_new / d * PERTURBATION_SIZE
end
-
+
return sum_log / (n_iterations * time_step)
end
export lyapunov_exponent
@@ -62,62 +67,68 @@ export lyapunov_exponent
Calculate the full Lyapunov spectrum for a continuous dynamical system using QR decomposition.
# Arguments
-- `model`: The dynamical model
-- `x0`: Initial condition vector
-- `time_step`: Total integration time per iteration
-- `dt`: Time step for integration (default: 0.01)
-- `warmup`: Number of warmup iterations (default: 1000)
-- `n_iterations`: Number of iterations for calculation (default: 10000)
+
+ - `model`: The dynamical model
+ - `x0`: Initial condition vector
+ - `time_step`: Total integration time per iteration
+ - `dt`: Time step for integration (default: 0.01)
+ - `warmup`: Number of warmup iterations (default: 1000)
+ - `n_iterations`: Number of iterations for calculation (default: 10000)
# Returns
-- Vector of Lyapunov exponents (sorted in descending order)
+
+ - Vector of Lyapunov exponents (sorted in descending order)
"""
-function lyapunov_spectrum(model, x0::AbstractVector{T}, time_step::Real;
- dt=0.01, warmup=1000, n_iterations=10000) where T
+function lyapunov_spectrum(
+ model, x0::AbstractVector{T}, time_step::Real;
+ dt = 0.01, warmup = 1000, n_iterations = 10000
+ ) where {T}
dim = length(x0)
x = copy(x0)
-
+
# Initialize orthonormal perturbation vectors
Q = Matrix{T}(I, dim, dim)
-
+
# Warmup
for _ in 1:warmup
x, _ = rk4_step(model, 0.0, x, time_step, dt)
# Evolve perturbation vectors
Q_evolved = zeros(T, dim, dim)
for j in 1:dim
- Q_evolved[:, j], _ = rk4_step(model, 0.0, x + Q[:, j] * PERTURBATION_SIZE, time_step, dt)
+ Q_evolved[:, j], _ =
+ rk4_step(model, 0.0, x + Q[:, j] * PERTURBATION_SIZE, time_step, dt)
Q_evolved[:, j] = Q_evolved[:, j] - x
end
Q, _ = qr(Q_evolved)
end
-
+
# Calculation
sum_log = zeros(T, dim)
for _ in 1:n_iterations
x_new, _ = rk4_step(model, 0.0, x, time_step, dt)
-
+
# Evolve perturbation vectors
Q_evolved = zeros(T, dim, dim)
for j in 1:dim
- Q_evolved[:, j], _ = rk4_step(model, 0.0, x + Q[:, j] * PERTURBATION_SIZE, time_step, dt)
+ Q_evolved[:, j], _ =
+ rk4_step(model, 0.0, x + Q[:, j] * PERTURBATION_SIZE, time_step, dt)
Q_evolved[:, j] = Q_evolved[:, j] - x_new
end
-
+
# QR decomposition
Q_new, R = qr(Q_evolved)
-
+
# Accumulate logarithms of diagonal elements
for j in 1:dim
sum_log[j] += log(abs(R[j, j]))
end
-
+
x = x_new
Q = Matrix(Q_new)
end
-
+
λs = sum_log / (n_iterations * time_step)
- return sort(λs, rev=true)
+ return sort(λs; rev = true)
end
export lyapunov_spectrum
@@ -128,20 +139,20 @@ Perform Runge-Kutta 4th order integration for time T with step dt.
Internal helper function for Lyapunov exponent calculations.
"""
-function rk4_step(f, t::Real, x::AbstractVector{T}, total_time::Real, dt::Real) where T
+function rk4_step(f, t::Real, x::AbstractVector{T}, total_time::Real, dt::Real) where {T}
n_steps = Int(round(total_time / dt))
current_t = t
current_x = copy(x)
-
+
for _ in 1:n_steps
k1 = f(current_t, current_x)
- k2 = f(current_t + dt/2, current_x + dt/2 * k1)
- k3 = f(current_t + dt/2, current_x + dt/2 * k2)
+ k2 = f(current_t + dt / 2, current_x + dt / 2 * k1)
+ k3 = f(current_t + dt / 2, current_x + dt / 2 * k2)
k4 = f(current_t + dt, current_x + dt * k3)
-
- current_x = current_x + (dt/6) * (k1 + 2*k2 + 2*k3 + k4)
+
+ current_x = current_x + (dt / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
current_t += dt
end
-
+
return current_x, current_t
end
diff --git a/src/analysis/poincare.jl b/src/analysis/poincare.jl
index fcee1ca..782450e 100644
--- a/src/analysis/poincare.jl
+++ b/src/analysis/poincare.jl
@@ -4,18 +4,21 @@
Calculate a Poincaré section for a continuous dynamical system.
# Arguments
-- `model`: The dynamical model (must be callable with signature `model(t, x)`)
-- `x0`: Initial condition vector
-- `plane_normal`: Normal vector to the Poincaré section plane
-- `plane_point`: A point on the Poincaré section plane
-- `t_max`: Maximum integration time
-- `dt`: Time step for integration (default: 0.01)
-- `direction`: Direction of crossing to record: `:positive`, `:negative`, or `:both` (default: `:both`)
+
+ - `model`: The dynamical model (must be callable with signature `model(t, x)`)
+ - `x0`: Initial condition vector
+ - `plane_normal`: Normal vector to the Poincaré section plane
+ - `plane_point`: A point on the Poincaré section plane
+ - `t_max`: Maximum integration time
+ - `dt`: Time step for integration (default: 0.01)
+ - `direction`: Direction of crossing to record: `:positive`, `:negative`, or `:both` (default: `:both`)
# Returns
-- Vector of points where trajectory crosses the Poincaré section
+
+ - Vector of points where trajectory crosses the Poincaré section
# Example
+
```julia
model = Lorenz()
x0 = [1.0, 1.0, 1.0]
@@ -25,42 +28,44 @@ plane_point = [0.0, 0.0, 27.0]
section = poincare_section(model, x0, plane_normal, plane_point, 1000.0)
```
"""
-function poincare_section(model, x0::AbstractVector{T},
- plane_normal::AbstractVector{T},
- plane_point::AbstractVector{T},
- t_max::Real;
- dt=0.01,
- direction=:both) where T
+function poincare_section(
+ model, x0::AbstractVector{T},
+ plane_normal::AbstractVector{T},
+ plane_point::AbstractVector{T},
+ t_max::Real;
+ dt = 0.01,
+ direction = :both
+ ) where {T}
# Normalize the plane normal
n = plane_normal / norm(plane_normal)
-
+
section_points = Vector{Vector{T}}()
-
+
t = 0.0
x = copy(x0)
-
+
# Calculate signed distance from plane
dist_prev = dot(x - plane_point, n)
-
+
n_steps = ceil(Int, t_max / dt)
-
+
for _ in 1:n_steps
# RK4 step
k1 = model(t, x)
- k2 = model(t + dt/2, x + dt/2 * k1)
- k3 = model(t + dt/2, x + dt/2 * k2)
+ k2 = model(t + dt / 2, x + dt / 2 * k1)
+ k3 = model(t + dt / 2, x + dt / 2 * k2)
k4 = model(t + dt, x + dt * k3)
-
- x_new = x + (dt/6) * (k1 + 2*k2 + 2*k3 + ...*[Comment body truncated]* |
Codecov Report❌ Patch coverage is
🚀 New features to boost your workflow:
|
- Add FormatFix.yml (replaces FormatCheck.yml): auto-commits Blue style formatting on PR - Add PRLabeler.yml: labels PR from checkbox selections in template - Add PublishRelease.yml: creates GitHub release on git tag push - Add VersionCheck.yml: enforces Project.toml version bump on every PR to main - Add AutoMerge.yml (updated): trigger on Bot PRs or automerge label - Add .github/release-drafter.yml and PULL_REQUEST_TEMPLATE.md - Update CI.yml, Documentation.yml: checkout@v6 and clean up comments Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Description
DataVault.jl という自作のモジュールに合わせてもろもろを変更しました
documentが中途半端だったりするので、その中途半端さをある程度解消していく予定です
Fixed: # (if any)
Type of Change
enhancement)bug)performance)documentation)choreorrefactor)Proposed Changes
what did you change?
Usage or Results
# Example or verification scriptcheck list