Home Python TutorialMatplotlib Library Python | Decay , Bayes Update ,Double Pendulum problem and Oscilloscope in Matplotlib

# Python | Decay , Bayes Update ,Double Pendulum problem and Oscilloscope in Matplotlib

## Decay

This example showcases : – using a generator to drive an animation, – changing axes limits during an animation.

``````import itertools
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def data_gen():
for cnt in itertools.count():
t = cnt / 10
yield t, np.sin(2*np.pi*t) * np.exp(-t/10)
def init():
ax.set_ylim(-1.1, 1.1)
ax.set_xlim(0, 10)
del xdata[:]
del ydata[:]
line.set_data(xdata, ydata)
return line,
fig, ax = plt.subplots()
line = ax.plot([], [], lw=2)
ax.grid()
xdata, ydata = [], []

def run(data):
# update the data
t, y = data
xdata.append(t)
ydata.append(y)
xmin, xmax = ax.get_xlim()

if t >= xmax:
ax.set_xlim(xmin, 2*xmax)
ax.figure.canyas.draw()
line.set_data(xdata, ydata)
return line,
ani = animation.FuncAnimation(fig, run, data_gen, interval=10,init_func=init)
plt.show()``````

## The Bayes update

This animation displays the posterior estimate updates as it is refitted when new data arrives. The vertical line represents the theoretical value to which the plotted distribution should converge.

``````import math

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def beta_pdf(x, a, b):
return (x**(a-1) * (1-x)**(b-1) * math.gamma(a + b)
/ (math.gamma(a) * math.gamma(b)))

class UpdateDist:
def __init__(self, ax, prob=0.5):
self.success = 0
self.prob = prob
self.line, = ax.plot([], [], 'k-')
self.x = np.linspace(0, 1, 200)
self.ax = ax

# Set up plot parameters
self.ax.set_xlim(0, 1)
self.ax.set_ylim(0, 10)
self.ax.grid(True)

# This vertical line represents the theoretical value, to
# which the plotted distribution should converge.
self.ax.axvline(prob, linestyle='--', color='black')

def __call__(self, i):
# This way the plot can continuously run and we just keep
# watching new realizations of the process
if i == 0:
self.success = 0
self.line.set_data([], [])
return self.line,

# Choose success based on exceed a threshold with a uniform pick
if np.random.rand(1,) < self.prob:
self.success += 1
y = beta_pdf(self.x, self.success + 1, (i - self.success) + 1)
self.line.set_data(self.x, y)
return self.line,

# Fixing random state for reproducibility
np.random.seed(19680801)

fig, ax = plt.subplots()
ud = UpdateDist(ax, prob=0.7)
anim = FuncAnimation(fig, ud, frames=100, interval=100, blit=True)
plt.show()``````

## The double pendulum problem

``````from numpy import sin, cos
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import matplotlib.animation as animation

G = 9.8  # acceleration due to gravity, in m/s^2
L1 = 1.0  # length of pendulum 1 in m
L2 = 1.0  # length of pendulum 2 in m
M1 = 1.0  # mass of pendulum 1 in kg
M2 = 1.0  # mass of pendulum 2 in kg
t_stop = 5  # how many seconds to simulate

def derivs(state, t):

dydx = np.zeros_like(state)
dydx[0] = state[1]

delta = state[2] - state[0]
den1 = (M1+M2) * L1 - M2 * L1 * cos(delta) * cos(delta)
dydx[1] = ((M2 * L1 * state[1] * state[1] * sin(delta) * cos(delta)
+ M2 * G * sin(state[2]) * cos(delta)
+ M2 * L2 * state[3] * state[3] * sin(delta)
- (M1+M2) * G * sin(state[0]))
/ den1)

dydx[2] = state[3]

den2 = (L2/L1) * den1
dydx[3] = ((- M2 * L2 * state[3] * state[3] * sin(delta) * cos(delta)
+ (M1+M2) * G * sin(state[0]) * cos(delta)
- (M1+M2) * L1 * state[1] * state[1] * sin(delta)
- (M1+M2) * G * sin(state[2]))
/ den2)

return dydx

# create a time array from 0..100 sampled at 0.05 second steps
dt = 0.05
t = np.arange(0, t_stop, dt)

# th1 and th2 are the initial angles (degrees)
# w10 and w20 are the initial angular velocities (degrees per second)
th1 = 120.0
w1 = 0.0
th2 = -10.0
w2 = 0.0

# initial state
state = np.radians([th1, w1, th2, w2])

# integrate your ODE using scipy.integrate.
y = integrate.odeint(derivs, state, t)

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

x2 = L2*sin(y[:, 2]) + x1
y2 = -L2*cos(y[:, 2]) + y1

fig = plt.figure(figsize=(5, 4))
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 1))
ax.set_aspect('equal')
ax.grid()

line, = ax.plot([], [], 'o-', lw=2)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

def animate(i):
thisx = [0, x1[i], x2[i]]
thisy = [0, y1[i], y2[i]]

line.set_data(thisx, thisy)
time_text.set_text(time_template % (i*dt))
return line, time_text

ani = animation.FuncAnimation(
fig, animate, len(y), interval=dt*1000, blit=True)
plt.show()``````

## Oscilloscope

``````import numpy as np
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt
import matplotlib.animation as animation

class Scope:
def __init__(self, ax, maxt=2, dt=0.02):
self.ax = ax
self.dt = dt
self.maxt = maxt
self.tdata = [0]
self.ydata = [0]
self.line = Line2D(self.tdata, self.ydata)
self.ax.set_ylim(-.1, 1.1)
self.ax.set_xlim(0, self.maxt)

def update(self, y):
lastt = self.tdata[-1]
if lastt > self.tdata[0] + self.maxt:  # reset the arrays
self.tdata = [self.tdata[-1]]
self.ydata = [self.ydata[-1]]
self.ax.set_xlim(self.tdata[0], self.tdata[0] + self.maxt)
self.ax.figure.canvas.draw()

t = self.tdata[-1] + self.dt
self.tdata.append(t)
self.ydata.append(y)
self.line.set_data(self.tdata, self.ydata)
return self.line,

def emitter(p=0.1):
"""Return a random value in [0, 1) with probability p, else 0."""
while True:
v = np.random.rand(1)
if v > p:
yield 0.
else:
yield np.random.rand(1)

# Fixing random state for reproducibility
np.random.seed(19680801 // 10)

fig, ax = plt.subplots()
scope = Scope(ax)

# pass a generator in "emitter" to produce data for the update func
ani = animation.FuncAnimation(fig, scope.update, emitter, interval=50,
blit=True)

plt.show()``````