-
Notifications
You must be signed in to change notification settings - Fork 64
Expand file tree
/
Copy pathphase_plot_ode.py
More file actions
137 lines (96 loc) · 2.89 KB
/
phase_plot_ode.py
File metadata and controls
137 lines (96 loc) · 2.89 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
"""
Phase portrait of ordinary differential equations
=================================================
Very basic example of how to plot the trajectories of a few trajectories that solve an ordinary differential equation
given some constants that help define the starting points
"""
# test_example = false
# sphinx_gallery_pygfx_docs = 'screenshot'
import numpy as np
import scipy
import fastplotlib as fpl
def create_and_solve_ode(
t_start: float,
t_end: float,
n_t: int,
p: float,
q: float,
A: np.ndarray,
coef: float = 1.0,
) -> float:
"""
Given an ODE of the form:
x' = Ax
where:
x' ∈ R^2
x ∈ R^2, x is a function of x, i.e. x(t)
A ∈ R^2x2
We can solve this ODE through an eigendecomposition of A:
A η_i = λ_i η_i
where λ_i, η_i are the eigenvalues and eigenvectors of A.
Thus this gives us the following solution for x(t):
x(t) = p * η_1 * e^(λ_1 * t) + q * η_2 * e^(λ_2 * t)
where p and q are constants that determine the initial conditions
References
----------
https://tutorial.math.lamar.edu/Classes/DE/SolutionsToSystems.aspx
Parameters
----------
t_start: float
t_end: float
p: float
q: float
A: np.ndarray
coef: float
Returns
-------
np.ndarray
"""
matrix = A * np.array([[1, 1], [1 / coef, 1 / coef]])
(lam1, lam2), evecs = scipy.linalg.eig(matrix)
eta1 = evecs[:, 0]
eta2 = evecs[:, 1]
x_t = np.zeros((n_t, 2))
for i, t in enumerate(np.linspace(t_start, t_end, n_t)):
x_t[i] = (p * eta1 * np.exp(lam1 * t)) + (q * eta2 * np.exp(lam2 * t))
return x_t
# constants range
START, STOP = -2, 2
T_START, T_END = -2, 5
A = np.array([
[0, 1],
[-1, -1]
])
x_t_collection = list()
x_t_vel_collection = list()
for p in np.linspace(START, STOP, 8):
for q in np.linspace(START, STOP, 8):
n_t = 100
x_t_vel = np.zeros(n_t)
x_t = create_and_solve_ode(
t_start=T_START,
t_end=T_END,
n_t=n_t,
p=p,
q=q,
A=A,
coef=0.5,
)
x_t_vel[:-1] = np.linalg.norm(x_t[:-1] - x_t[1:], axis=1, ord=2)
x_t_vel[-1] = x_t_vel[-2]
x_t_collection.append(x_t)
x_t_vel_collection.append(x_t_vel)
fig = fpl.Figure()
lines = fig[0, 0].add_line_collection(x_t_collection, cmap="viridis", cmap_transform=x_t_vel_collection)
for (line, vel) in zip(lines, x_t_vel_collection):
line.cmap = "viridis"
# color indicates velocity along trajectory
line.cmap.transform = vel
print(vel)
fig[0, 0].axes.intersection = (0, 0, 0)
fig.show()
# NOTE: `if __name__ == "__main__"` is NOT how to use fastplotlib interactively
# please see our docs for using fastplotlib interactively in ipython and jupyter
if __name__ == "__main__":
print(__doc__)
fpl.loop.run()