|
1 | 1 | # Example: Defining a `RateSystem` |
2 | 2 |
|
| 3 | +## Basic concept |
| 4 | + |
3 | 5 | Consider an autonomous deterministic dynamical system `ds` (e.g. a `CoupledODEs`) of which |
4 | 6 | you want to ramp one parameter, i.e. change the parameter's value over time. |
5 | 7 |
|
6 | | -Using the `RateSystem` type, you can easily achieve this in a two-step process: |
7 | | -1) Specify a forcing profile `p(t)` over an interval ``t\in I``. This profile, stored as a `ForcingProfile` type, |
8 | | - describes the shape of the parameter ramping you would like to consider. |
9 | | -2) Apply this parametric forcing to a given autonomous dynamical system by constructing a `RateSystem`. |
10 | | - This yields a non-autonomous dynamical system in which the parameters are explicitly time-dependent. |
| 8 | +Using the `RateSystem` type, you can easily achieve this in two steps: |
| 9 | +1. Specify a [`ForcingProfile`](@ref) that describes the shape of the parameter ramping `p(t)` over an interval ``t\in I`` |
| 10 | +2. Apply this parametric forcing to the system `ds` by constructing a [`RateSystem`](@ref), i.e. a nonautonomous system in which the parameters are explicitly time-dependent. |
| 11 | + |
| 12 | + |
11 | 13 |
|
12 | 14 | You can rescale the forcing profile in system time units by specifying the start time and |
13 | | -duration of the forcing change. Then: |
14 | | -- for times `t < forcing_start_time`, the system is autonomous, with parameters given by the underlying autonomous system |
15 | | -- for `forcing_start_time < t < forcing_start_time + forcing_duration`, the system is non-autonomous with the parameter change given by the `ForcingProfile` |
16 | | -- for `t > forcing_start_time + forcing_duration`, the system is again autonomous, with parameters fixed at their values attained at the end of the forcing interval (i.e. `t = forcing_start_time + forcing_duration`). |
| 15 | +duration of the forcing change. Then, for |
| 16 | +- `t < forcing_start_time` |
| 17 | + - the system is autonomous, with parameters given by the underlying autonomous system |
| 18 | +- `forcing_start_time < t < forcing_start_time + forcing_duration` |
| 19 | + - the system is non-autonomous with the parameter change given by the `ForcingProfile`, scaled in magnitude by `forcing_scale` |
| 20 | +- `t > forcing_start_time + forcing_duration` |
| 21 | + - the system is again autonomous, with parameters fixed at their values attained at the end of the forcing interval (i.e. `t = forcing_start_time + forcing_duration`). |
| 22 | + |
17 | 23 | This setting is widely used and convenient for studying rate-dependent tipping. |
18 | 24 |
|
19 | 25 | ## Example |
20 | 26 |
|
21 | | -Let us explore a simple prototypical example. |
22 | | - |
23 | | -We consider the following one-dimensional autonomous system with one attractor, given by |
| 27 | +As a simple prototypical example, let's consider the following one-dimensional autonomous system with one stable fixed point, given by |
24 | 28 | the ordinary differential equation: |
25 | 29 | ```math |
26 | 30 | \begin{aligned} |
27 | 31 | \dot{x} &= (x+p)^2 - 1 |
28 | 32 | \end{aligned} |
29 | 33 | ``` |
30 | | -The parameter ``p`` shifts the location of the extrema of the drift field. |
| 34 | +The parameter ``p`` shifts the location of the equilibria ``\dot x = 0``. |
31 | 35 | We implement this system as follows: |
32 | 36 |
|
33 | | -````julia |
| 37 | +```@example MAIN |
34 | 38 | using CriticalTransitions |
35 | | -using CairoMakie |
36 | 39 |
|
37 | | -function f(u, p, t) # out-of-place |
| 40 | +function f(u, p, t) |
38 | 41 | x = u[1] |
39 | 42 | dx = (x + p[1])^2 - 1 |
40 | 43 | return SVector{1}(dx) |
41 | 44 | end |
42 | 45 |
|
43 | | -x0 = [-1.0] |
44 | | -ds = CoupledODEs(f, x0, [0.0]) |
45 | | -```` |
| 46 | +x0 = [-1.0] # Initial state |
| 47 | +p0 = [0.0] # Initial parameter value |
46 | 48 |
|
47 | | -## Applying the parameter ramping |
| 49 | +ds = CoupledODEs(f, x0, p0) # Autonomous system |
| 50 | +``` |
48 | 51 |
|
49 | 52 | Now, we want to explore a non-autonomous version of this system by applying a parameter |
50 | | -shift s.t. speed and amplitude of the parameter shift are easily modifiable but the 'shape' |
51 | | -of it is always the same - just linearly stretched or squeezed. |
| 53 | +change over a given time interval. |
| 54 | + |
| 55 | +### Forcing profile |
| 56 | +First, create a `ForcingProfile` to specify the functional form of the parameter change `p(t)`, here a hyperbolic tangent: |
| 57 | + |
| 58 | +```@example MAIN |
| 59 | +profile(t) = tanh(t) |
| 60 | +interval = (-5.0, 5.0) |
| 61 | +
|
| 62 | +fp = ForcingProfile(profile, interval) |
| 63 | +``` |
52 | 64 |
|
53 | | -First specify a section of a function `p(t)` that you would like to use to ramp a |
54 | | -parameter of `ds`: |
| 65 | +Let's plot the forcing profile: |
55 | 66 |
|
56 | | -````julia |
57 | | -p(t) = tanh(t) # A monotonic function that describes the parameter shift |
58 | | -interval = (-5.0, 5.0) # Domain interval of p(t) we want to consider |
59 | | -fp = ForcingProfile(p, interval) |
60 | | -```` |
| 67 | +```@example MAIN |
| 68 | +using CairoMakie |
| 69 | +
|
| 70 | +time_range = range(fp.interval[1], fp.interval[2]; length=100); |
| 71 | +
|
| 72 | +fig = Figure(); |
| 73 | +ax = Axis(fig[1, 1]; xlabel="t (arbitrary units)", ylabel="profile (arbitrary units)"); |
| 74 | +lines!(ax, time_range, fp.profile.(time_range); linewidth=2); |
| 75 | +fig |
| 76 | +``` |
61 | 77 |
|
62 | | -Then specify how you would like to use the `ForcingProfile` to shift the `pidx`'th parameter |
63 | | -of ds: |
| 78 | +Note that the `interval` is given in arbitrary units - the profile is rescaled to your system's units in the next step. |
64 | 79 |
|
65 | | -````julia |
66 | | -pidx = 1 # Index of the parameter within the parameter-container of ds |
67 | | -forcing_start_time = -50.0 # Time when the parameter shift should start |
68 | | -forcing_duration = 105.0 # Time interval over which p(interval) is spread out or squeezed |
69 | | -forcing_scale = 3.0 # Amplitude of the ramping. `p` is then automatically rescaled |
70 | | -t0 = -70.0 # Initial time of the resulting non-autonomous system (relevant to later compute trajectories) |
| 80 | +### Applying the forcing |
71 | 81 |
|
72 | | -rs = RateSystem(ds, fp, pidx; forcing_start_time, forcing_duration, forcing_scale, t0) |
| 82 | +Now, specify how the forcing profile `fp` should be applied to the `pidx`-th parameter of your system `ds` by constructing a `RateSystem`. |
73 | 83 |
|
74 | | -#note # Choosing different values of the `forcing_duration` allows us to vary the speed of the parameter ramping, while its shape remains the same, and it only gets stretched or squeezed. |
| 84 | +```@example MAIN |
| 85 | +pidx = 1 # parameter index |
| 86 | +forcing_start_time = 20.0 # system time units |
| 87 | +forcing_duration = 105.0 # system time units |
| 88 | +forcing_scale = 3.0 |
| 89 | +t0 = 0.0 # system initial time |
75 | 90 |
|
76 | | -#note # If `p(t)` within the `ForcingProfile` is a monotonic function, the `forcing_scale` will give the total amplitude of the parameter ramping. For non-monotonic `p(t)`, the `forcing_scale` will only linearly scale the amplitude of the parameter ramping, but does not equal the total amplitude. |
77 | | -```` |
| 91 | +rs = RateSystem(ds, fp, pidx; |
| 92 | + forcing_start_time, forcing_duration, forcing_scale, t0) |
| 93 | +``` |
78 | 94 |
|
79 | | -Now, we can compute trajectories of this new system `rate_system` and of the previous autonomous system `ds` in the familiar way: |
| 95 | +The `forcing_scale` is a multiplication factor that scales the profile `fp.profile`. Here, we have ``p(5)-p(-5) \approx 2``, so the amplitude of the parameter change is ``6`` after multiplying with `forcing_scale = 3`. |
80 | 96 |
|
81 | | -````julia |
82 | | -T = forcing_duration + 40.0; # length of the trajectory that we want to compute |
83 | | -auto_traj = trajectory(ds, T, x0); |
84 | | -nonauto_traj = trajectory(rs, T, x0); |
85 | | -```` |
| 97 | +In the `RateSystem`, the time dependence of the parameter `p[pidx]` thus looks like this: |
86 | 98 |
|
87 | | -We plot the two trajectories: |
| 99 | +```@example MAIN |
| 100 | +T = forcing_duration + 40.0 # Total time |
| 101 | +t_points = range(t0, t0+T, length=100) |
| 102 | +
|
| 103 | +parameter(t) = parameters(rs, t)[pidx] # Returns the parameter value at time t |
88 | 104 |
|
89 | | -````julia |
90 | 105 | fig = Figure(); |
91 | | -axs = Axis(fig[1, 1]; xlabel="t", ylabel="x"); |
| 106 | +ax = Axis(fig[1, 1]; xlabel="Time (system units)", ylabel="p[1]"); |
| 107 | +lines!(ax, t_points, parameter.(t_points); linewidth=2); |
| 108 | +fig |
| 109 | +``` |
| 110 | + |
| 111 | +!!! tip "Modifying the forcing" |
| 112 | + |
| 113 | + In a `RateSystem`, the forcing can easily be modified to implement different forcing rates and forcing amplitudes. |
| 114 | + - Via `set_forcing_duration!(rs, length)`, you can change the length of the forcing interval and thus the rate of change of the forcing. |
| 115 | + - Via `set_forcing_scale!(rs, factor)`, you can change the magnitude of the forcing by a given factor. |
| 116 | + |
| 117 | +### Simulating trajectories |
| 118 | + |
| 119 | +The `RateSystem` type behaves just like the type of underlyinh autonomous system, in this case a `CoupledODEs`. Thus, we can simply call the `trajectory` function to simulate either the autonomous system `ds` or the nonautonomous system `rs`. |
| 120 | + |
| 121 | +```@example MAIN |
| 122 | +traj_ds = trajectory(ds, T, x0) |
| 123 | +traj_rs = trajectory(rs, T, x0) |
| 124 | +``` |
| 125 | + |
| 126 | +Let's compare the two trajectories: |
| 127 | + |
| 128 | +```@example MAIN |
| 129 | +fig = Figure(); |
| 130 | +axs = Axis(fig[1, 1]; xlabel="Time", ylabel="x"); |
92 | 131 | lines!( |
93 | 132 | axs, |
94 | | - t0 .+ auto_traj[2], |
95 | | - auto_traj[1][:, 1]; |
| 133 | + t0 .+ traj_ds[2], |
| 134 | + traj_ds[1][:, 1]; |
96 | 135 | linewidth=2, |
97 | | - label=L"\text{Autonomous system}", |
| 136 | + label="Autonomous CoupledODEs (ds)", |
98 | 137 | ); |
99 | 138 | lines!( |
100 | 139 | axs, |
101 | | - nonauto_traj[2], |
102 | | - nonauto_traj[1][:, 1]; |
| 140 | + traj_rs[2], |
| 141 | + traj_rs[1][:, 1]; |
103 | 142 | linewidth=2, |
104 | | - label=L"\text{Nonautonomous system}", |
| 143 | + label="Nonautonomous RateSystem (rs)", |
105 | 144 | ); |
106 | | -axislegend(axs; position=:rc, labelsize=10); |
| 145 | +axislegend(axs; position=:lb); |
107 | 146 | fig |
108 | | -```` |
109 | | - |
110 | | -We can also plot the shifted parameter `p(t)`: |
111 | | - |
112 | | -````julia |
113 | | -time_range = range(t0, t0+T; length=Int(T+1)); |
114 | | - |
115 | | -fig = Figure(); |
116 | | -ax = Axis(fig[1, 1]; xlabel="t", ylabel="p(t)"); |
117 | | -lines!(ax, time_range, fp.profile.(time_range); linewidth=2); |
118 | | -fig |
119 | | -```` |
120 | | - |
121 | | ---- |
122 | | - |
123 | | -*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* |
| 147 | +``` |
124 | 148 |
|
| 149 | +While the autonomous system `ds` remains at the fixed point ``x^*=-1``, the nonautonomous system tracks the moving equilibrium until reaching the stable fixed point ``x^*=-7`` of the future limit system (i.e. the autonomous limit system after the parameter change) where ``p=6``. |
0 commit comments