Skip to content

Commit

Permalink
style: rename x0 to be u0 in line with SciML (#290)
Browse files Browse the repository at this point in the history
  • Loading branch information
oameye authored Oct 26, 2024
1 parent 26cf081 commit d91edcb
Show file tree
Hide file tree
Showing 4 changed files with 19 additions and 19 deletions.
14 changes: 7 additions & 7 deletions docs/src/tutorials/time_dependent.md
Original file line number Diff line number Diff line change
Expand Up @@ -43,10 +43,10 @@ Given $\mathbf{u}(T_0)$, what is $\mathbf{u}(T)$ at future times?
For constant parameters, a [`HarmonicEquation`](@ref HarmonicBalance.HarmonicEquation) object can be fed into the constructor of [`ODEProblem`](@ref ODEProblem). The syntax is similar to DifferentialEquations.jl :
```@example time_dependent
using OrdinaryDiffEqTsit5
x0 = [0.; 0.] # initial condition
u0 = [0.; 0.] # initial condition
fixed = (ω0 => 1.0, γ => 1e-2, λ => 5e-2, F => 1e-3, α => 1.0, η => 0.3, θ => 0, ω => 1.0) # parameter values
ode_problem = ODEProblem(harmonic_eq, fixed, x0 = x0, timespan = (0,1000))
ode_problem = ODEProblem(harmonic_eq, fixed, u0 = u0, timespan = (0,1000))
```
OrdinaryDiffEq.jl takes it from here - we only need to use `solve`.

Expand All @@ -55,10 +55,10 @@ time_evo = solve(ode_problem, Tsit5(), saveat=1.0);
plot(time_evo, ["u1", "v1"], harmonic_eq)
```

Running the above code with `x0 = [0.2, 0.2]` gives the plots
Running the above code with `u0 = [0.2, 0.2]` gives the plots
```@example time_dependent
x0 = [0.2; 0.2] # initial condition
ode_problem = remake(ode_problem, u0 = x0)
u0 = [0.2; 0.2] # initial condition
ode_problem = remake(ode_problem, u0 = u0)
time_evo = solve(ode_problem, Tsit5(), saveat=1.0);
plot(time_evo, ["u1", "v1"], harmonic_eq)
```
Expand All @@ -70,7 +70,7 @@ result = get_steady_states(harmonic_eq, varied, fixed)
plot(result, "sqrt(u1^2 + v1^2)")
```

Clearly when evolving from `x0 = [0.,0.]`, the system ends up in the low-amplitude branch 2. With `x0 = [0.2, 0.2]`, the system ends up in branch 3.
Clearly when evolving from `u0 = [0., 0.]`, the system ends up in the low-amplitude branch 2. With `u0 = [0.2, 0.2]`, the system ends up in branch 3.

## Adiabatic parameter sweeps

Expand All @@ -84,7 +84,7 @@ The sweep linearly interpolates between $\omega = 0.9$ at time 0 and $\omega =

Let us now define a new `ODEProblem` which incorporates `sweep` and again use `solve`:
```@example time_dependent
ode_problem = ODEProblem(harmonic_eq, fixed, sweep=sweep, x0=[0.1;0.0], timespan=(0, 2e4))
ode_problem = ODEProblem(harmonic_eq, fixed, sweep=sweep, u0=[0.1;0.0], timespan=(0, 2e4))
time_evo = solve(ode_problem, Tsit5(), saveat=100)
plot(time_evo, "sqrt(u1^2 + v1^2)", harmonic_eq)
```
Expand Down
16 changes: 8 additions & 8 deletions ext/TimeEvolution/ODEProblem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,20 +4,20 @@ using OrdinaryDiffEqTsit5: ODEProblem, solve, ODESolution
ODEProblem(
eom::HarmonicEquation;
fixed_parameters,
x0::Vector,
sweep::AdiabaticSweep,
u0::Vector,
sweep::ParameterSweep,
timespan::Tuple
)
Creates an ODEProblem object used by OrdinaryDiffEqTsit5.jl from the equations in `eom` to simulate time-evolution within `timespan`.
`fixed_parameters` must be a dictionary mapping parameters+variables to numbers (possible to use a solution index, e.g. solutions[x][y] for branch y of solution x).
If `x0` is specified, it is used as an initial condition; otherwise the values from `fixed_parameters` are used.
If `u0` is specified, it is used as an initial condition; otherwise the values from `fixed_parameters` are used.
"""
function OrdinaryDiffEqTsit5.ODEProblem(
eom::HarmonicEquation,
fixed_parameters;
sweep::AdiabaticSweep=AdiabaticSweep(),
x0::Vector=[],
sweep::ParameterSweep=ParameterSweep(),
u0::Vector=[],
timespan::Tuple,
perturb_initial=0.0,
kwargs...,
Expand Down Expand Up @@ -54,11 +54,11 @@ function OrdinaryDiffEqTsit5.ODEProblem(
end
end

# the initial condition is x0 if specified, taken from fixed_parameters otherwise
initial = if isempty(x0)
# the initial condition is u0 if specified, taken from fixed_parameters otherwise
initial = if isempty(u0)
real.(collect(values(fixed_parameters))[1:length(vars)]) * (1 - perturb_initial)
else
x0
u0
end

return ODEProblem(f!, initial, timespan; kwargs...)
Expand Down
4 changes: 2 additions & 2 deletions test/SteadyStateDiffEqExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,8 @@ using HarmonicBalance,
model = structural_simplify(model)

param = [α, ω, ω0, F, γ] .=> [1.0, 1.2, 1.0, 0.01, 0.01]
x0 = [1.0, 0.0]
prob_ss = SteadyStateProblem{true}(model, x0, param; jac=true)
u0 = [1.0, 0.0]
prob_ss = SteadyStateProblem{true}(model, u0, param; jac=true)
prob_np = NonlinearProblem(prob_ss)

ω_span = (0.9, 1.5)
Expand Down
4 changes: 2 additions & 2 deletions test/time_evolution.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,8 @@ varied = ω => range(0.9, 1.1, 10)
res = get_steady_states(p, varied, fixed; show_progress=false, seed=SEED)

tspan = (0.0, 10)
sweep = AdiabaticSweep=> (0.9, 1.1), tspan) # linearly interpolate between two values at two times
ode_problem = ODEProblem(harmonic_eq, fixed; sweep=sweep, x0=[0.01; 0.0], timespan=tspan)
sweep = ParameterSweep=> (0.9, 1.1), tspan) # linearly interpolate between two values at two times
ode_problem = ODEProblem(harmonic_eq, fixed; sweep=sweep, u0=[0.01; 0.0], timespan=tspan)
time_soln = solve(ode_problem, Tsit5(); saveat=1);

transform_solutions(time_soln, "sqrt(u1^2+v1^2)", harmonic_eq)

0 comments on commit d91edcb

Please sign in to comment.