# The second order differential equation for the angle `theta` of a
# pendulum acted on by gravity with friction can be written::

# theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0

# where `b` and `c` are positive constants, and a prime (') denotes a
# derivative.  To solve this equation with `odeint`, we must first convert
# it to a system of first order equations.  By defining the angular
# velocity ``omega(t) = theta'(t)``, we obtain the system::

# theta'(t) = omega(t)
# omega'(t) = -b*omega(t) - c*sin(theta(t))

# Let `y` be the vector [`theta`, `omega`].  We implement this system
# in python as:

def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt
# ...

# We assume the constants are `b` = 0.25 and `c` = 5.0:

b = 0.25
c = 5.0

# For initial conditions, we assume the pendulum is nearly vertical
# with `theta(0)` = `pi` - 0.1, and it initially at rest, so
# `omega(0)` = 0.  Then the vector of initial conditions is

y0 = [np.pi - 0.1, 0.0]

# We generate a solution 101 evenly spaced samples in the interval
# 0 <= `t` <= 10.  So our array of times is:

t = np.linspace(0, 10, 101)

# Call `odeint` to generate the solution.  To pass the parameters
# `b` and `c` to `pend`, we give them to `odeint` using the `args`
# argument.

from scipy.integrate import odeint
sol = odeint(pend, y0, t, args=(b, c))

# The solution is an array with shape (101, 2).  The first column
# is `theta(t)`, and the second is `omega(t)`.  The following code
# plots both components.

import matplotlib.pyplot as plt
plt.plot(t, sol[:, 0], 'b', label='theta(t)')
plt.plot(t, sol[:, 1], 'g', label='omega(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()