Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

Appearance settings

Commit e67514f

Browse filesBrowse files
committed
Cleanup differential equations examples.
Slightly cleanup the lorenz_attractor example; also use Euler's method for the double pendulum. It is less accurate, but sufficient for illustrative purposes (also, 1. the double pendulum is chaotic anyways so even RK4 will end up quite far from the actual behavior, if one waits for long enough, as can be checked by trying various integrators in solve_ivp; 2. we're fine with using Euler's method for the also chaotic lorenz_attractor). The point is also to make fewer examples dependent on scipy.
1 parent a35921c commit e67514f
Copy full SHA for e67514f

File tree

2 files changed

+31
-29
lines changed
Filter options

2 files changed

+31
-29
lines changed

‎examples/animation/double_pendulum.py

Copy file name to clipboardExpand all lines: examples/animation/double_pendulum.py
+13-7Lines changed: 13 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,6 @@
1212
from numpy import sin, cos
1313
import numpy as np
1414
import matplotlib.pyplot as plt
15-
import scipy.integrate as integrate
1615
import matplotlib.animation as animation
1716
from collections import deque
1817

@@ -22,13 +21,13 @@
2221
L = L1 + L2 # maximal length of the combined pendulum
2322
M1 = 1.0 # mass of pendulum 1 in kg
2423
M2 = 1.0 # mass of pendulum 2 in kg
25-
t_stop = 5 # how many seconds to simulate
24+
t_stop = 2.5 # how many seconds to simulate
2625
history_len = 500 # how many trajectory points to display
2726

2827

29-
def derivs(state, t):
30-
28+
def derivs(t, state):
3129
dydx = np.zeros_like(state)
30+
3231
dydx[0] = state[1]
3332

3433
delta = state[2] - state[0]
@@ -51,7 +50,7 @@ def derivs(state, t):
5150
return dydx
5251

5352
# create a time array from 0..t_stop sampled at 0.02 second steps
54-
dt = 0.02
53+
dt = 0.01
5554
t = np.arange(0, t_stop, dt)
5655

5756
# th1 and th2 are the initial angles (degrees)
@@ -64,8 +63,15 @@ def derivs(state, t):
6463
# initial state
6564
state = np.radians([th1, w1, th2, w2])
6665

67-
# integrate your ODE using scipy.integrate.
68-
y = integrate.odeint(derivs, state, t)
66+
# integrate the ODE using Euler's method
67+
y = np.empty((len(t), 4))
68+
y[0] = state
69+
for i in range(1, len(t)):
70+
y[i] = y[i - 1] + derivs(t[i - 1], y[i - 1]) * dt
71+
72+
# A more accurate estimate could be obtained e.g. using scipy:
73+
#
74+
# y = scipy.integrate.solve_ivp(derivs, t[[0, -1]], state, t_eval=t).y.T
6975

7076
x1 = L1*sin(y[:, 0])
7177
y1 = -L1*cos(y[:, 0])

‎examples/mplot3d/lorenz_attractor.py

Copy file name to clipboardExpand all lines: examples/mplot3d/lorenz_attractor.py
+18-22Lines changed: 18 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -18,45 +18,41 @@
1818
import matplotlib.pyplot as plt
1919

2020

21-
def lorenz(x, y, z, s=10, r=28, b=2.667):
21+
def lorenz(xyz, *, s=10, r=28, b=2.667):
2222
"""
23-
Given:
24-
x, y, z: a point of interest in three dimensional space
25-
s, r, b: parameters defining the lorenz attractor
26-
Returns:
27-
x_dot, y_dot, z_dot: values of the lorenz attractor's partial
28-
derivatives at the point x, y, z
23+
Parameters
24+
----------
25+
xyz : array-like, shape (3,)
26+
Point of interest in three dimensional space.
27+
s, r, b : float
28+
Parameters defining the Lorenz attractor.
29+
30+
Returns
31+
-------
32+
xyz_dot : array, shape (3,)
33+
Values of the Lorenz attractor's partial derivatives at *xyz*.
2934
"""
35+
x, y, z = xyz
3036
x_dot = s*(y - x)
3137
y_dot = r*x - y - x*z
3238
z_dot = x*y - b*z
33-
return x_dot, y_dot, z_dot
39+
return np.array([x_dot, y_dot, z_dot])
3440

3541

3642
dt = 0.01
3743
num_steps = 10000
3844

39-
# Need one more for the initial values
40-
xs = np.empty(num_steps + 1)
41-
ys = np.empty(num_steps + 1)
42-
zs = np.empty(num_steps + 1)
43-
44-
# Set initial values
45-
xs[0], ys[0], zs[0] = (0., 1., 1.05)
46-
45+
xyzs = np.empty((num_steps + 1, 3)) # Need one more for the initial values
46+
xyzs[0] = (0., 1., 1.05) # Set initial values
4747
# Step through "time", calculating the partial derivatives at the current point
4848
# and using them to estimate the next point
4949
for i in range(num_steps):
50-
x_dot, y_dot, z_dot = lorenz(xs[i], ys[i], zs[i])
51-
xs[i + 1] = xs[i] + (x_dot * dt)
52-
ys[i + 1] = ys[i] + (y_dot * dt)
53-
zs[i + 1] = zs[i] + (z_dot * dt)
54-
50+
xyzs[i + 1] = xyzs[i] + lorenz(xyzs[i]) * dt
5551

5652
# Plot
5753
ax = plt.figure().add_subplot(projection='3d')
5854

59-
ax.plot(xs, ys, zs, lw=0.5)
55+
ax.plot(*xyzs.T, lw=0.5)
6056
ax.set_xlabel("X Axis")
6157
ax.set_ylabel("Y Axis")
6258
ax.set_zlabel("Z Axis")

0 commit comments

Comments
0 (0)
Morty Proxy This is a proxified and sanitized view of the page, visit original site.