text:
# Day 13

```python
! cat README.md
```

```python
with open('input.txt', 'rt') as f:
    timestamp = int(f.readline().rstrip())
    services = [int(x) for x in f.readline().split(',') if x != 'x']
```

```python
def partone(timestamp, services):
    mindiff = min([(s - (timestamp % s), s) for s in services], key=lambda x: x[0])
    return mindiff[0]*mindiff[1]
```

```python
# test
test_timestamp = 939
test_services = [7,13,59,31,19]
partone(test_timestamp, test_services)
```

```python
# solution
partone(timestamp, services)
```
--- Part Two ---
The shuttle company is running a contest: one gold coin for anyone that can find the earliest timestamp such that the first bus ID departs at that time and each subsequent listed bus ID departs at that subsequent minute. (The first line in your input is no longer relevant.)

For example, suppose you have the same list of bus IDs as above:

7,13,x,x,59,x,31,19
An x in the schedule means there are no constraints on what bus IDs must depart at that time.

This means you are looking for the earliest timestamp (called t) such that:

Bus ID 7 departs at timestamp t.
Bus ID 13 departs one minute after timestamp t.
There are no requirements or restrictions on departures at two or three minutes after timestamp t.
Bus ID 59 departs four minutes after timestamp t.
There are no requirements or restrictions on departures at five minutes after timestamp t.
Bus ID 31 departs six minutes after timestamp t.
Bus ID 19 departs seven minutes after timestamp t.
The only bus departures that matter are the listed bus IDs at their specific offsets from t. Those bus IDs can depart at other times, and other bus IDs can depart at those times. For example, in the list above, because bus ID 19 must depart seven minutes after the timestamp at which bus ID 7 departs, bus ID 7 will always also be departing with bus ID 19 at seven minutes after timestamp t.

In this example, the earliest timestamp at which this occurs is 1068781:

time     bus 7   bus 13  bus 59  bus 31  bus 19
1068773    .       .       .       .       .
1068774    D       .       .       .       .
1068775    .       .       .       .       .
1068776    .       .       .       .       .
1068777    .       .       .       .       .
1068778    .       .       .       .       .
1068779    .       .       .       .       .
1068780    .       .       .       .       .
1068781    D       .       .       .       .
1068782    .       D       .       .       .
1068783    .       .       .       .       .
1068784    .       .       .       .       .
1068785    .       .       D       .       .
1068786    .       .       .       .       .
1068787    .       .       .       D       .
1068788    D       .       .       .       D
1068789    .       .       .       .       .
1068790    .       .       .       .       .
1068791    .       .       .       .       .
1068792    .       .       .       .       .
1068793    .       .       .       .       .
1068794    .       .       .       .       .
1068795    D       D       .       .       .
1068796    .       .       .       .       .
1068797    .       .       .       .       .
In the above example, bus ID 7 departs at timestamp 1068788 (seven minutes after t). This is fine; the only requirement on that minute is that bus ID 19 departs then, and it does.

Here are some other examples:

The earliest timestamp that matches the list 17,x,13,19 is 3417.
67,7,59,61 first occurs at timestamp 754018.
67,x,7,59,61 first occurs at timestamp 779210.
67,7,x,59,61 first occurs at timestamp 1261476.
1789,37,47,1889 first occurs at timestamp 1202161486.
However, with so many bus IDs in your list, surely the actual earliest timestamp will be larger than 100000000000000!

What is the earliest timestamp such that all of the listed bus IDs depart at offsets matching their positions in the list?

```python
with open('input.txt', 'rt') as f:
    timestamp = int(f.readline().rstrip())
    services = [x for x in f.readline().rstrip().split(',')]
```

```python
from sympy.ntheory.modular import solve_congruence

def parttwo(services):
    rem_mod_pairs = [((-1) * i, int(s)) for i, s in enumerate(services) if s != 'x']
    return int(solve_congruence(*rem_mod_pairs)[0])
```

```python
# test
test_services = ['7','13','x','x','59','x','31','19']
parttwo(test_services)
```

```python
# solution
parttwo(services)
```

max_stars_repo_name: bbglab/adventofcode

max_issues_repo_name: bbglab/adventofcode

max_forks_repo_name: bbglab/adventofcode


text:
```python
from tectosaur2.nb_config import setup

setup()
```

# Quasidynamic earthquake cycles on two faults

## Summary

* Much like the last section, we simulate quasidynamic rupture evolution but this time on two faults, one of which is curved.
* Most of the code here follows almost exactly from [the previous section on strike-slip/antiplane earthquake cycles](part6_qd).
* Realistic tectonic loading is a challenge when there is more than one fault.
* The boundary integral approach used here is barely more complex with two faults than with one fault.

## Loading a multi-fault system

In the single fault system from the last section, we were able to gloss over loading by using a slip deficit or backslip formulation. Because the physical system is invariant with respect to a rigid body displacement, any tectonic motion that has already been fully accomodated by fault slip can be ignored. Put differently, in a single fault system with a purely elastic rheology, the stresses from tectonic motion can be equivalently described by backwards slip on the fault.

\begin{equation}
u_{\textrm{total}} = u_{\textrm{loading}} + u_{\textrm{fault}}
\end{equation}

There are a number of different options for loading a fault system:
- **backslip**: as discussed, we simply apply a backwards slip rate to account for distant plate motion. Once the stresses from this backwards slip accumulate, the fault slips forward.
- **farfield displacement**: apply a constant velocity boundary condition at a distance from the fault.
- **stress**: Either on a lower boundary of the domain or on distant side boundaries, apply a constant traction/stress loading rate.
- **block model**: this is a more advanced version of backslip loading where each fault-bounded subregion of the domain has its own block rotation/motion rate. These rates are assumed to be the long term slip rate. A corresponding kinematic slip deficit is applied to each fault.

Here, I do something a bit like the block model solution by just applying a backslip of 31 mm of loading per year (1 nm of loading per second) to both faults equally. This would model a real-world situation with a plate boundary that has 62 mm per year of strike-slip motion. Exactly 50% of that motion is accomodated by each of the two faults.

```python
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from tectosaur2 import (
    gauss_rule,
    integrate_term,
    panelize_symbolic_surface,
    concat_meshes,
)
from tectosaur2.laplace2d import double_layer, hypersingular

```

```python
shear_modulus = 3.2e10

quad_rule = gauss_rule(5)
t = sp.var("t")

fault_bottom = 40000
fault1 = panelize_symbolic_surface(
    t, -5000 + t * 0, fault_bottom * (t + 1) * -0.5,
    quad_rule,
    n_panels=200
)

fault2_x = 5000 + (t * 10000 + 10000) + (t ** 2 * 5000 - 5000)
fault2 = panelize_symbolic_surface(
    t, fault2_x, fault_bottom * (t + 1) * -0.5,
    quad_rule,
    n_panels=200
)
faults = concat_meshes((fault1, fault2))

surf_half_L = 100000
free = panelize_symbolic_surface(
    t, -t * surf_half_L, 0 * t,
    quad_rule,
    n_panels=200
)
```

```python
print(
    f"The free surface mesh has {free.n_panels} panels with a total of {free.n_pts} points."
)
print(
    f"The fault mesh has {faults.n_panels} panels with a total of {faults.n_pts} points."
)
```

    The free surface mesh has 200 panels with a total of 1000 points.
    The fault mesh has 400 panels with a total of 2000 points.

```python
plt.plot(free.pts[:,0]/1000, free.pts[:,1]/1000, 'k-')
plt.plot(fault1.pts[:,0]/1000, fault1.pts[:,1]/1000, 'r-')
plt.plot(fault2.pts[:,0]/1000, fault2.pts[:,1]/1000, 'r-')
plt.xlabel(r'$x ~ \mathrm{(km)}$')
plt.ylabel(r'$y ~ \mathrm{(km)}$')
plt.axis('scaled')
plt.xlim([-100, 100])
plt.ylim([-80, 20])
plt.show()
```

```python
singularities = np.array(
    [
        [-surf_half_L, 0],
        [surf_half_L, 0],
        [5000, 0],
        [-5000, 0]
    ]
)

(free_disp_to_free_disp, fault_slip_to_free_disp), report = integrate_term(
    double_layer, free.pts, free, faults, tol=1e-11, singularities=singularities, return_report=True
)
fault_slip_to_free_disp = fault_slip_to_free_disp[:, 0, :, 0]
free_disp_solve_mat = (
    np.eye(free_disp_to_free_disp.shape[0]) + free_disp_to_free_disp[:, 0, :, 0]
)
```

```python
(free_disp_to_fault_stress, fault_slip_to_fault_stress), report = integrate_term(
    hypersingular,
    faults.pts,
    free,
    faults,
    tol=1e-10,
    singularities=singularities,
    return_report=True,
)
fault_slip_to_fault_stress *= shear_modulus
free_disp_to_fault_stress *= shear_modulus

fault_slip_to_fault_traction = np.sum(
    fault_slip_to_fault_stress[:,:,:,0] * faults.normals[:, :, None], axis=1
)
free_disp_to_fault_traction = np.sum(
    free_disp_to_fault_stress[:,:,:,0] * faults.normals[:, :, None], axis=1
)
```

    /Users/tbent/Dropbox/active/eq/tectosaur2/tectosaur2/integrate.py:208: UserWarning: Some expanded integrals reached maximum expansion order. These integrals may be inaccurate.
      warnings.warn(

```python
A = -fault_slip_to_fault_traction
B = -free_disp_to_fault_traction
C = fault_slip_to_free_disp
Dinv = np.linalg.inv(free_disp_solve_mat)

total_fault_slip_to_fault_traction = A - B.dot(Dinv.dot(C))
```

```python
analytical_sd = -np.arctan(-fault_bottom / free.pts[:,0]) / np.pi
numerical_sd = Dinv.dot(-C).sum(axis=1)
rp = (faults.pts[:,1] + 1) ** 2
analytical_ft = -(shear_modulus / (2 * np.pi)) * ((1.0 / (faults.pts[:,1] + fault_bottom)) - (1.0 / (faults.pts[:,1] - fault_bottom)))
numerical_ft = total_fault_slip_to_fault_traction.sum(axis=1)
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.title('Stress')
plt.plot(faults.pts[:,1], np.log10(np.abs(numerical_ft)), 'b-')
plt.xlim([-fault_bottom, 0])
plt.xlabel('$\mathrm{depth ~~ (m)}$')
plt.ylabel('$t$')
plt.subplot(1,2,2)
plt.title('Displacement')
plt.plot(free.pts[:,0], numerical_sd, 'b-')
plt.xlim([-100000, 100000])
plt.xlabel('$x ~~ \mathrm{(m)}$')
plt.ylabel('$u$')
plt.show()
```

## Rate and state friction

```python
siay = 31556952 # seconds in a year
density = 2670  # rock density (kg/m^3)
cs = np.sqrt(shear_modulus / density)  # Shear wave speed (m/s)
eta = shear_modulus / (2 * cs)  # The radiation damping coefficient (kg / (m^2 * s))
Vp = 1e-9  # Rate of plate motion
sigma_n = 50e6  # Normal stress (Pa)

# parameters describing "a", the coefficient of the direct velocity strengthening effect
a0 = 0.01
amax = 0.025
fy = faults.pts[:, 1]
H = 15000
h = 3000
fp_a = np.where(
    fy > -H, a0, np.where(fy > -(H + h), a0 + (amax - a0) * (fy + H) / -h, amax)
)

# state-based velocity weakening effect
fp_b = 0.015

# state evolution length scale (m)
fp_Dc = 0.008

# baseline coefficient of friction
fp_f0 = 0.6

# if V = V0 and at steady state, then f = f0, units are m/s
fp_V0 = 1e-6
```

```python
plt.figure(figsize=(3, 5))
plt.plot(fp_a, fy/1000, label='a')
plt.plot(np.full(fy.shape[0], fp_b), fy/1000, label='b')
plt.xlim([0, 0.03])
plt.ylabel('depth')
plt.legend()
plt.show()
```

```python
# The regularized form of the aging law.
def aging_law(fp_a, V, state):
    return (fp_b * fp_V0 / fp_Dc) * (np.exp((fp_f0 - state) / fp_b) - (V / fp_V0))

def qd_equation(fp_a, shear_stress, V, state):
    # The regularized rate and state friction equation
    F = sigma_n * fp_a * np.arcsinh(V / (2 * fp_V0) * np.exp(state / fp_a))

    # The full shear stress balance:
    return shear_stress - eta * V - F

def qd_equation_dV(fp_a, V, state):
    # First calculate the derivative of the friction law with respect to velocity
    # This is helpful for equation solving using Newton's method
    expsa = np.exp(state / fp_a)
    Q = (V * expsa) / (2 * fp_V0)
    dFdV = fp_a * expsa * sigma_n / (2 * fp_V0 * np.sqrt(1 + Q * Q))

    # The derivative of the full shear stress balance.
    return -eta - dFdV

def solve_friction(fp_a, shear, V_old, state, tol=1e-10):
    V = V_old
    max_iter = 150
    for i in range(max_iter):
        # Newton's method step!
        f = qd_equation(fp_a, shear, V, state)
        dfdv = qd_equation_dV(fp_a, V, state)
        step = f / dfdv

        # We know that slip velocity should not be negative so any step that
        # would lead to negative slip velocities is "bad". In those cases, we
        # cut the step size in half iteratively until the step is no longer
        # bad. This is a backtracking line search.
        while True:
            bad = step > V
            if np.any(bad):
                step[bad] *= 0.5
            else:
                break

        # Take the step and check if everything has converged.
        Vn = V - step
        if np.max(np.abs(step) / Vn) < tol:
            break
        V = Vn
        if i == max_iter - 1:
            raise Exception("Failed to converge.")

    # Return the solution and the number of iterations required.
    return Vn, i
```

```python
for fm in [fault1, fault2]:
    mesh_L = np.max(np.abs(np.diff(fm.pts[:, 1])))
    Lb = shear_modulus * fp_Dc / (sigma_n * fp_b)
    hstar = (np.pi * shear_modulus * fp_Dc) / (sigma_n * (fp_b - fp_a))
    print(mesh_L, Lb, np.min(hstar[hstar > 0]))
```

    53.84693101057201 341.3333333333333 3216.990877275949
    53.84693101057201 341.3333333333333 3216.990877275949

## Quasidynamic earthquake cycle derivatives

```python
from scipy.optimize import fsolve
import copy

init_state_scalar = fsolve(lambda S: aging_law(fp_a, Vp, S), 0.7)[0]
init_state = np.full(faults.n_pts, init_state_scalar)
tau_amax = -qd_equation(amax, 0, Vp, init_state_scalar)

init_traction = np.full(faults.n_pts, tau_amax)
init_slip_deficit = np.zeros(faults.n_pts)
init_conditions = np.concatenate((init_slip_deficit, init_state))
```

```python
class SystemState:

    V_old = np.full(faults.n_pts, Vp)
    state = None

    def calc(self, t, y, verbose=False):
        # Separate the slip_deficit and state sub components of the
        # time integration state.
        slip_deficit = y[: init_slip_deficit.shape[0]]
        state = y[init_slip_deficit.shape[0] :]

        # If the state values are bad, then the adaptive integrator probably
        # took a bad step.
        if np.any((state < 0) | (state > 2.0)):
            return False

        # The big three lines solving for quasistatic shear stress, slip rate
        # and state evolution
        tau_qs = init_traction - total_fault_slip_to_fault_traction.dot(slip_deficit)
        V = solve_friction(fp_a, tau_qs, self.V_old, state)[0]
        dstatedt = aging_law(fp_a, V, state)
        if not np.all(np.isfinite(V)):
            return False
        self.V_old = V

        slip_deficit_rate = Vp - V
        out = slip_deficit, state, tau_qs, V, slip_deficit_rate, dstatedt
        self.state = out
        return self.state

def plot_system_state(SS):
    """This is just a helper function that creates some rough plots of the
    current state to help with debugging"""
    slip_deficit, state, tau_qs, slip_deficit_rate, dstatedt = SS

    plt.figure(figsize=(15, 9))
    plt.subplot(2, 3, 1)
    plt.title("slip")
    plt.plot(faults.pts[:, 1], slip_deficit)

    plt.subplot(2, 3, 2)
    plt.title("state")
    plt.plot(faults.pts[:, 1], state)

    plt.subplot(2, 3, 3)
    plt.title("tau qs")
    plt.plot(faults.pts[:, 1], tau_qs)

    plt.subplot(2, 3, 4)
    plt.title("slip rate")
    plt.plot(faults.pts[:, 1], slip_deficit_rate)

    plt.subplot(2, 3, 6)
    plt.title("dstatedt")
    plt.plot(faults.pts[:, 1], dstatedt)
    plt.tight_layout()

    plt.show()

def calc_derivatives(state, t, y):
    """
    This helper function calculates the system state and then extracts the
    relevant derivatives that the integrator needs. It also intentionally
    returns infinite derivatives when the `y` vector provided by the integrator
    is invalid.
    """
    if not np.all(np.isfinite(y)):
        return np.inf * y
    state = state.calc(t, y)
    if not state:
        return np.inf * y
    derivatives = np.concatenate((state[-2], state[-1]))
    return derivatives
```

## Integrating through time

```python
%%time
from scipy.integrate import RK45

# We use a 5th order adaptive Runge Kutta method and pass the derivative function to it
# the relative tolerance will be 1e-11 to make sure that even
state = SystemState()
derivs = lambda t, y: calc_derivatives(state, t, y)
integrator = RK45
atol = Vp * 1e-6
rtol = 1e-11
rk = integrator(derivs, 0, init_conditions, 1e50, atol=atol, rtol=rtol)

# Set the initial time step to one day.
rk.h_abs = 60 * 60 * 24

# Integrate for 1000 years.
max_T = 1000 * siay

n_steps = 1000000
t_history = [0]
y_history = [init_conditions.copy()]
for i in range(n_steps):

    # Take a time step and store the result
    if rk.step() != None:
        raise Exception("TIME STEPPING FAILED")
        break
    t_history.append(rk.t)
    y_history.append(rk.y.copy())

    # Print the time every 5000 steps
    if i % 5000 == 0:
        print(f"step={i}, time={rk.t / siay} yrs, step={(rk.t - t_history[-2]) / siay}")
    if rk.t > max_T:
        break

y_history = np.array(y_history)
t_history = np.array(t_history)
```

    step=0, time=0.002737907006988508 yrs, step=0.002737907006988508
    step=5000, time=161.2388012289919 yrs, step=6.068323795478569e-11
    step=10000, time=161.238801483585 yrs, step=3.200375945922213e-11
    step=15000, time=177.06735503165402 yrs, step=0.05782401363757997
    step=20000, time=187.60702659014385 yrs, step=8.29256996752649e-11
    step=25000, time=187.6070269264095 yrs, step=2.786049654906221e-10
    step=30000, time=213.4153146695212 yrs, step=1.3348499105890855e-10
    step=35000, time=229.49165911338014 yrs, step=1.8458825569115088e-10
    step=40000, time=298.5155537603046 yrs, step=1.1030568652139828e-10
    step=45000, time=298.51555399934733 yrs, step=3.251751197178755e-11
    step=50000, time=298.5168884817506 yrs, step=6.51359998844683e-06
    step=55000, time=323.01280183052285 yrs, step=1.2281707123916784e-10
    step=60000, time=323.01280207887106 yrs, step=3.3484575524851865e-11
    step=65000, time=323.2817940109216 yrs, step=0.001729718400292269
    step=70000, time=351.52626665410605 yrs, step=1.0785478482910092e-08
    step=75000, time=368.9516542294202 yrs, step=1.2976784052681758e-10
    step=80000, time=368.951654876993 yrs, step=6.515590688770813e-11
    step=85000, time=428.9915445603336 yrs, step=0.0014206634043130257
    step=90000, time=446.2939275656298 yrs, step=6.588120455250637e-11
    step=95000, time=446.29392801977554 yrs, step=9.456068304806992e-10
    step=100000, time=469.8994523682489 yrs, step=1.0256917809688377e-10
    step=105000, time=469.89945268723017 yrs, step=5.1798341561007304e-11
    step=110000, time=469.9812354270467 yrs, step=0.00034767164051762437
    step=115000, time=500.90614938849876 yrs, step=1.439715864624497e-10
    step=120000, time=573.9064649517662 yrs, step=1.424001081887202e-10
    step=125000, time=573.9064651906277 yrs, step=2.5022769435539118e-11
    step=130000, time=597.5674917048766 yrs, step=0.0012973735857483035
    step=135000, time=597.6230358626846 yrs, step=5.391379308333549e-11
    step=140000, time=597.6230366893858 yrs, step=4.587749495737104e-09
    step=145000, time=625.2112327038584 yrs, step=1.1906969997104363e-10
    step=150000, time=642.3455004578122 yrs, step=1.2898210138995284e-10
    step=155000, time=704.8417989867936 yrs, step=0.003113930973116129
    step=160000, time=716.3612963754118 yrs, step=2.9978970144993726e-11
    step=165000, time=716.3612966987445 yrs, step=3.9843018386249724e-10
    step=170000, time=743.0069768079197 yrs, step=1.6222491102653862e-10
    step=175000, time=743.0069772376432 yrs, step=4.6660816435353134e-11
    step=180000, time=743.0070823319467 yrs, step=6.836717438689674e-07
    step=185000, time=772.4229461837367 yrs, step=1.1205848921132735e-10
    step=190000, time=790.5550252030326 yrs, step=9.598105764163313e-11
    step=195000, time=790.5550257931098 yrs, step=9.694812119469745e-11
    step=200000, time=790.5770082962302 yrs, step=0.00011880289499414454
    step=205000, time=843.4093589520743 yrs, step=1.3986156636192637e-10
    step=210000, time=843.4093595399053 yrs, step=5.7540281407326665e-11
    step=215000, time=865.8601485323139 yrs, step=4.5113514750450234e-10
    step=220000, time=903.1263077569464 yrs, step=1.2438854951289735e-10
    step=225000, time=903.1263082887026 yrs, step=8.413452911659529e-11
    step=230000, time=927.4523165504868 yrs, step=0.008045238903279458
    step=235000, time=947.0642174585033 yrs, step=1.7528026899290684e-11
    step=240000, time=947.0685830652857 yrs, step=1.768826485736451e-05
    step=245000, time=972.5636479272499 yrs, step=4.078590535048743e-10
    step=250000, time=988.09794803525 yrs, step=8.195863612220059e-11
    step=255000, time=988.0979484325295 yrs, step=2.7077779485800782e-11
    CPU times: user 2h 15min 3s, sys: 36min 52s, total: 2h 51min 55s
    Wall time: 1h 48min 42s

## Plotting the results

```python
derivs_history = np.diff(y_history, axis=0) / np.diff(t_history)[:, None]
max_vel = np.max(np.abs(derivs_history), axis=1)
plt.plot(t_history[1:] / siay, np.log10(max_vel))
plt.xlabel('$t ~~ \mathrm{(yrs)}$')
plt.ylabel('$\log_{10}(V)$')
plt.show()
```

```python
for j in range(2):
    fm = fault1 if j == 0 else fault2
    plt.figure(figsize=(10, 4))
    last_plt_t = -1000
    last_plt_slip = init_slip_deficit
    event_times = []
    for i in range(len(y_history) - 1):
        y = y_history[i]
        t = t_history[i]
        slip_deficit_full = y[: init_slip_deficit.shape[0]]
        fault1_sd = slip_deficit_full[:fault1.n_pts]
        fault2_sd = slip_deficit_full[fault1.n_pts:]
        slip_deficit = fault1_sd if j == 0 else fault2_sd
        should_plot = False

        # Plot a red line every three second if the slip rate is over 0.1 mm/s.
        if (
            max_vel[i] >= 0.0001 and t - last_plt_t > 3
        ):
            if len(event_times) == 0 or t - event_times[-1] > siay:
                event_times.append(t)
            should_plot = True
            color = "r"

        # Plot a blue line every fifteen years during the interseismic period
        if t - last_plt_t > 15 * siay:
            should_plot = True
            color = "b"

        if should_plot:
            # Convert from slip deficit to slip:
            slip = -slip_deficit + Vp * t
            plt.plot(slip, fm.pts[:,1] / 1000.0, color + "-", linewidth=0.5)
            last_plt_t = t
            last_plt_slip = slip
    plt.xlim([0, np.max(last_plt_slip)])
    plt.ylim([-40, 0])
    plt.ylabel(r"$\textrm{z (km)}$")
    plt.xlabel(r"$\textrm{slip (m)}$")
    plt.tight_layout()
    plt.savefig("halfspace.png", dpi=300)
    plt.show()
```

ext: ipynb


text:
```python
%matplotlib inline
```

Autograd: Automatic Differentiation
===================================

Central to all neural networks in PyTorch is the ``autograd`` package.
Let’s first briefly visit this, and we will then go to training our
first neural network.

The ``autograd`` package provides automatic differentiation for all operations
on Tensors. It is a define-by-run framework, which means that your backprop is
defined by how your code is run, and that every single iteration can be
different.

Let us see this in more simple terms with some examples.

Tensor
--------

``torch.Tensor`` is the central class of the package. If you set its attribute
``.requires_grad`` as ``True``, it starts to track all operations on it. When
you finish your computation you can call ``.backward()`` and have all the
gradients computed automatically. The gradient for this tensor will be
accumulated into ``.grad`` attribute.

To stop a tensor from tracking history, you can call ``.detach()`` to detach
it from the computation history, and to prevent future computation from being
tracked.

To prevent tracking history (and using memory), you can also wrap the code block
in ``with torch.no_grad():``. This can be particularly helpful when evaluating a
model because the model may have trainable parameters with
``requires_grad=True``, but for which we don't need the gradients.

There’s one more class which is very important for autograd
implementation - a ``Function``.

``Tensor`` and ``Function`` are interconnected and build up an acyclic
graph, that encodes a complete history of computation. Each tensor has
a ``.grad_fn`` attribute that references a ``Function`` that has created
the ``Tensor`` (except for Tensors created by the user - their
``grad_fn is None``).

If you want to compute the derivatives, you can call ``.backward()`` on
a ``Tensor``. If ``Tensor`` is a scalar (i.e. it holds a one element
data), you don’t need to specify any arguments to ``backward()``,
however if it has more elements, you need to specify a ``gradient``
argument that is a tensor of matching shape.

```python
import torch
```

Create a tensor and set ``requires_grad=True`` to track computation with it

```python
x = torch.ones(2, 2, requires_grad=True)
print(x)
```

Do a tensor operation:

```python
y = x + 2
print(y)
```

``y`` was created as a result of an operation, so it has a ``grad_fn``.

```python
print(y.grad_fn)
```

Do more operations on ``y``

```python
z = y * y * 3
out = z.mean()

print(z, out)
```

``.requires_grad_( ... )`` changes an existing Tensor's ``requires_grad``
flag in-place. The input flag defaults to ``False`` if not given.

```python
a = torch.randn(2, 2)
a = ((a * 3) / (a - 1))
print(a.requires_grad)
a.requires_grad_(True)
print(a.requires_grad)
b = (a * a).sum()
print(b.grad_fn)
```

Gradients
---------
Let's backprop now.
Because ``out`` contains a single scalar, ``out.backward()`` is
equivalent to ``out.backward(torch.tensor(1.))``.

```python
out.backward()
```

Print gradients d(out)/dx

```python
print(x.grad)
```

You should have got a matrix of ``4.5``. Let’s call the ``out``
*Tensor* “$o$”.
We have that $o = \frac{1}{4}\sum_i z_i$,
$z_i = 3(x_i+2)^2$ and $z_i\bigr\rvert_{x_i=1} = 27$.
Therefore,
$\frac{\partial o}{\partial x_i} = \frac{3}{2}(x_i+2)$, hence
$\frac{\partial o}{\partial x_i}\bigr\rvert_{x_i=1} = \frac{9}{2} = 4.5$.

Mathematically, if you have a vector valued function $\vec{y}=f(\vec{x})$,
then the gradient of $\vec{y}$ with respect to $\vec{x}$
is a Jacobian matrix:

\begin{align}J=\left(\begin{array}{ccc}
   \frac{\partial y_{1}}{\partial x_{1}} & \cdots & \frac{\partial y_{1}}{\partial x_{n}}\\
   \vdots & \ddots & \vdots\\
   \frac{\partial y_{m}}{\partial x_{1}} & \cdots & \frac{\partial y_{m}}{\partial x_{n}}
   \end{array}\right)\end{align}

Generally speaking, ``torch.autograd`` is an engine for computing
vector-Jacobian product. That is, given any vector
$v=\left(\begin{array}{cccc} v_{1} & v_{2} & \cdots & v_{m}\end{array}\right)^{T}$,
compute the product $v^{T}\cdot J$. If $v$ happens to be
the gradient of a scalar function $l=g\left(\vec{y}\right)$,
that is,
$v=\left(\begin{array}{ccc}\frac{\partial l}{\partial y_{1}} & \cdots & \frac{\partial l}{\partial y_{m}}\end{array}\right)^{T}$,
then by the chain rule, the vector-Jacobian product would be the
gradient of $l$ with respect to $\vec{x}$:

\begin{align}J^{T}\cdot v=\left(\begin{array}{ccc}
   \frac{\partial y_{1}}{\partial x_{1}} & \cdots & \frac{\partial y_{m}}{\partial x_{1}}\\
   \vdots & \ddots & \vdots\\
   \frac{\partial y_{1}}{\partial x_{n}} & \cdots & \frac{\partial y_{m}}{\partial x_{n}}
   \end{array}\right)\left(\begin{array}{c}
   \frac{\partial l}{\partial y_{1}}\\
   \vdots\\
   \frac{\partial l}{\partial y_{m}}
   \end{array}\right)=\left(\begin{array}{c}
   \frac{\partial l}{\partial x_{1}}\\
   \vdots\\
   \frac{\partial l}{\partial x_{n}}
   \end{array}\right)\end{align}

(Note that $v^{T}\cdot J$ gives a row vector which can be
treated as a column vector by taking $J^{T}\cdot v$.)

This characteristic of vector-Jacobian product makes it very
convenient to feed external gradients into a model that has
non-scalar output.

Now let's take a look at an example of vector-Jacobian product:

```python
x = torch.randn(3, requires_grad=True)

y = x * 2
while y.data.norm() < 1000:
    y = y * 2

print(y)
```

Now in this case ``y`` is no longer a scalar. ``torch.autograd``
could not compute the full Jacobian directly, but if we just
want the vector-Jacobian product, simply pass the vector to
``backward`` as argument:

```python
v = torch.tensor([0.1, 1.0, 0.0001], dtype=torch.float)
y.backward(v)

print(x.grad)
```

You can also stop autograd from tracking history on Tensors
with ``.requires_grad=True`` by wrapping the code block in
``with torch.no_grad():``

```python
print(x.requires_grad)
print((x ** 2).requires_grad)

with torch.no_grad():
	print((x ** 2).requires_grad)
```

**Read Later:**

Document about ``autograd.Function`` is at
https://pytorch.org/docs/stable/autograd.html#function


text:
# OOP in Action: The Samuelson Accelerator

**Co-author: Natasha Watkins**

## Paul Samuelson’s (1939) multiplier-accelerator model

This lecture creates nonstochastic and stochastic versions of Paul
Samuelson’s celebrated multiplier accelerator model, as described in
[[Sam39]](zreferences.ipynb#samuelson1939)

Samuelson used a **second-order linear difference equation** to
represent a model of national output based on three components:

- a **national output identity** asserting that national outcome is the
  sum of output devoted to consumption plus investment plus government
  purchases
- a Keynesian **consumption function** asserting that consumption at
  time $ t $ is equal to a constant called the **marginal
  propensity to consume** times national output at time $ t-1 $
- an investment **accelerator** asserting that investment at time
  $ t $ equals a constant called the **accelerator coefficient**
  times the difference in output between period $ t-1 $ and
  $ t-2 $
- the idea that consumption plus investment plus government purchases
  constitute **aggregate demand** that automatically call forth an
  equal amount of **aggregate supply**

(To read about linear difference equations see
[here](https://en.wikipedia.org/wiki/Linear_difference_equation) or
chapter IX of [[Sar87]](zreferences.ipynb#sargent1987))

Samuelson used the model to analyze how particular values of the
marginal propensity to consume and the accelerator coefficient might
give rise to transient **business cycles** in national output with
alternative dynamic properties:

- smooth convergence to a constant level of output
- damped business cycles that eventually converge to a constant level
  of output
- persistent business cycles that neither dampen nor explode

Later in the lecture we present an extension to Samuelson’s model that
simply adds a random shock to the right side of the national income
identity (the random shock represents random fluctuations in aggregate demand

This modification makes national output become governed by a second-order
**stochastic linear difference equation** that, with appropriate parameter values,
gives rise to recurrent irregular business cycles.

(To read about stochastic linear difference equations see chapter XI of
[[Sar87]](zreferences.ipynb#sargent1987))

### Structure of the model

In more detail, let’s assume that

- $ \{G_t\} $ is a sequence of levels of government expenditures.
  We’ll start by setting $ G_t = G \forall t $
- $ \{C_t\} $ is a sequence of levels of aggregate consumption
  expenditures, a key endogenous variable in the model
- $ \{I_t\} $ is a sequence of rates of investment, another key
  endogenous variable
- $ \{Y_t\} $ is a sequence of levels of national income, yet
  another endogenous variable
- $ a $ is the marginal propensity to consume in the Keynesian
  consumption function $ C_t = a Y_{t-1} + \gamma $
- $ b $ is the “accelerator coefficient” in the “investment
  accerator” $I_t = b (Y_{t-1} - Y_{t-2}) $
- $ \{\epsilon_{t}\} $ is an IID sequence standard normal random variables
- $ \sigma \geq 0 $ is a “volatility”
  parameter — setting $ \sigma = 0 $ recovers the nonstochastic case
  that we’ll start with

The model combines the consumption function

$$
C_t = a Y_{t-1} + \gamma
$$

	(1)

with the investment accelerator

$$
I_t = b (Y_{t-1} - Y_{t-2})
$$

	(2)

and the national income identity

$$
Y_t = C_t + I_t + G_t
$$

	(3)

- The parameter $ a $ is peoples’ *marginal propensity to consume*
  out of income - equation [(1)](#equation-consumption) asserts that people consume a fraction of
  math:a in (0,1) of each additional dollar of income
- The parameter $b > 0 $ is the investment accelerator coefficient - equation
  [(2)](#equation-accelerator) asserts that people invest in physical capital when
  income is increasing and disinvest when it is decreasing

Equations [(1)](#equation-consumption), [(2)](#equation-accelerator), and [(3)](#equation-income_identity)
imply the following second-order linear difference equation for national income:

$$
Y_t = (a+b) Y_{t-1} - b Y_{t-2} + (\gamma + G_t)
$$

or

$$
Y_t = \rho_1 Y_{t-1} + \rho_2 Y_{t-2} + (\gamma + G_t)
$$

	(4)

where $ \rho_1 = (a+b) $ and $ \rho_2 = -b $

To complete the model, we require two **initial conditions**

If the model is to generate time series for $ t=0, \ldots, T $, we
require initial values

$$
Y_{-1} = \bar Y_{-1}, \quad  Y_{-2} = \bar Y_{-2}
$$

We’ll ordinarily set the parameters $ (a,b) $ so that starting from
an arbitrary pair of initial conditions
$ (\bar Y_{-1}, \bar Y_{-2}) $, national income $Y_t $ converges to
a constant value as $ t $ becomes large

We are interested in studying

- the transient fluctuations in $ Y_t $ as it converges to its
  **steady state** level
- the **rate** at which it converges to a steady state level

The deterministic version of the model described so far — meaning that
no random shocks hit aggregate demand — has only transient fluctuations

We can convert the model to one that has persistent irregular
fluctuations by adding a random shock to aggregate demand

### Stochastic version of the model

We create a **random** or **stochastic** version of the model by adding
a random process of **shocks** or **disturbances**
$ \{\sigma \epsilon_t \} $ to the right side of equation [(4)](#equation-second_order),
leading to the **second-order scalar linear stochastic difference
equation**:

$$
Y_t = G_t + a (1-b) Y_{t-1} - a b Y_{t-2} + \sigma \epsilon_{t}
$$

	(5)

### Mathematical analysis of the model

To get started, let’s set $ G_t \equiv 0 $, $ \sigma = 0 $, and
$ \gamma = 0 $.

Then we can write equation [(5)](#equation-second_stochastic) as

$$
Y_t = \rho_1 Y_{t-1} + \rho_2 Y_{t-2}
$$

or

$$
Y_{t+2} - \rho_1 Y_{t+1} - \rho_2 Y_t  = 0
$$

	(6)

To discover the properties of the solution of [(6)](#equation-second_stochastic2),
it is useful first to form the **characteristic polynomial**
for [(6)](#equation-second_stochastic2):

$$
z^2 - \rho_1 z  - \rho_2
$$

	(7)

where $ z $ is possibly a complex number

We want to find the two **zeros** (a.k.a. **roots**) – namely
$ \lambda_1, \lambda_2 $ – of the characteristic polynomial

These are two special values of $ z $, say $ z= \lambda_1 $ and
$ z= \lambda_2 $, such that if we set $ z $ equal to one of
these values in expresssion [(7)](#equation-polynomial),
the characteristic polynomial [(7)](#equation-polynomial) equals zero:

$$
z^2 - \rho_1 z  - \rho_2  = (z- \lambda_1 ) (z -\lambda_2) = 0
$$

	(8)

Equation [(8)](#equation-polynomial_sol) is said to **factor** the characteristic polynomial

When the roots are complex, they will occur as a complex conjugate pair

When the roots are complex, it is convenient to represent them in the
polar form

$$
\lambda_1 = r e^{i \theta}, \  \lambda_2 = r e^{-i \theta}
$$

where $ r $ is the *amplitude* of the complex number and
$ \theta $ is its *angle*.

These can also be represented as

$$
\lambda_1 = r (cos (\theta) + i \sin (\theta))
$$

$$
\lambda_2 = r (cos (\theta) - i \sin(\theta))
$$

(To read about the polar form, see
[here](https://www.varsitytutors.com/hotmath/hotmath_help/topics/polar-form-of-a-complex-number))

Given **initial conditions** $ Y_{-1}, Y_{-2} $ we want to generate
a **solution** of the difference equation [(6)](#equation-second_stochastic2)

It can be represented as

$$
Y_t = \lambda_1^t c_1 + \lambda_2^t c_2
$$

where $ c_1 $ and $ c_2 $ are constants that depend on the two
initial conditions and on $ \rho_1, \rho_2 $

When the roots are complex, some algebra that exploits the fact that the
roots appear as a complex conjugate pair implies that

$$
Y_t = \tilde c_1 r^t \cos(\theta t + \tilde c_1)
$$

where $ \tilde c_1, \tilde c_2 $ is a pair of constants chosen to
satisfy the given initial conditions for $ Y_{-1}, Y_{-2} $

This formula shows that when the roots are complex, $ Y_t $ displays
oscillations with **period** $ \check p =
\frac{2 \pi}{\theta} $ and **damping factor** $ r $

**Remark:** Following [[Sam39]](zreferences.ipynb#samuelson1939), we want to choose the parameters
$ a, b $ of the model so that the absolute values (of the possibly
complex) roots $ \lambda_1, \lambda_2 $ of the characteristic
polynomial are both strictly less than one:

$$
| \lambda_j | < 1 \quad \quad \text{for } j = 1, 2
$$

**Remark:** When both eigenvalues $ \lambda_1, \lambda_2 $ have
absolute values strictly less than one, the absolute value of the larger
one governs the rate of convergence to the steady state of the non
stochastic version of the model

### Things this lecture does

1. Writes a function to generate simulations of a $ \{Y_t\} $
   sequence as a function of time. Say about 40 or 80 periods. The
   function checks that $ a, b $ are set so that
   $ \lambda\_1, \lambda\_2 $ are less than
   unity in absolute value (also called “modulus”). The function also
   tells us whether the roots are complex, and, if they are complex,
   returns both their real and complex parts; if the roots are both
   real, the function returns their values. The function requires that
   we put in initial conditions for $ Y_{-1}, Y_{-2} $
1. The lecture produces a version of a graph in chapter IX, p. 189, of
   [[Sar87]](zreferences.ipynb#sargent1987). The
   graph maps values of $ \rho_1, \rho_2 $ into various regions –
   including ones that produce complex roots and therefore “cycles” in
   the outputs. We can use versions of this graph to help us guide how
   to set $ a,b $
1. We use our function written to simulate paths that are stochastic
   (when $ \sigma >0 $)
1. We have written the function in a way that allows us to input
   $ \{G_t\} $ paths of a few simple forms – e.g., one time jumps
   in $ G $ at some time; and a permanent jump in $ G $ that
   occurs at some time
1. We proceed to use the Samuelson multiplier-accererator model as a
   laboratory to make a simple OOP example that eventually can follow
   the example of the Solow model class in [this lecture](https://lectures.quantecon.org/py/python_oop.html).
   One lesson from
   doing this is that in writing the function and class to simulate the
   model, the “state” that determines next period’s $ Y_{t+1} $ is
   now not just the current value $ Y_t $ but also the once lagged
   value $ Y_{t-1} $. This involves a little more bookkeeping than is
   required in the Solow model class definition
1. We use the Samuelson multiplier-accelerator model as a vehicle for
   teaching how we can gradually add more and more fun features to the
   class. We want to have a method in the class that automatically
   generates a simulation, either nonstochastic ($ \sigma=0 $) or
   stochastic ($ \sigma > 0 $)
1. We also show how to map the Samuelson model into a simple instance of
   the LinearStateSpace class described [here](https://lectures.quantecon.org/py/linear_models.html).
   We can use a LinearStateSpace instance to do various things that we did above with
   our homemade function and class. Among other things, we show by
   example that the eigenvalues of the matrix $ A $ that we use to
   form the instance of the LinearStateSpace class for the Samuelson
   model equal the roots of the characteristic polynomial [(7)](#equation-polynomial) for the
   Samuelson multiplier accelerator model. Here is the formula for the
   matrix $ A $ in the linear state space system in the case that
   government expenditures are a constant $ G $:

$$
A = \begin{bmatrix} 1 & 0 & 0 \cr
                      \gamma + G & \rho_1 & \rho_2 \cr
                      0 & 1 & 0 \end{bmatrix}
$$

## Let’s get to work

We’ll start by drawing an informative graph from page 189 of [[Sar87]](zreferences.ipynb#sargent1987)

```python
import numpy as np
import matplotlib.pyplot as plt

def param_plot():

    """this function creates the graph on page 189 of Sargent Macroeconomic Theory, second edition, 1987"""

    fig, ax = plt.subplots(figsize=(12, 8))
    ax.set_aspect('equal')

    # Set axis
    xmin, ymin = -3, -2
    xmax, ymax = -xmin, -ymin
    plt.axis([xmin, xmax, ymin, ymax])

    # Set axis labels
    ax.set(xticks=[], yticks=[])
    ax.set_xlabel(r'$\rho_2$', fontsize=16)
    ax.xaxis.set_label_position('top')
    ax.set_ylabel(r'$\rho_1$', rotation=0, fontsize=16)
    ax.yaxis.set_label_position('right')

    # Draw (t1, t2) points
    rho1 = np.linspace(-2, 2, 100)
    ax.plot(rho1, -abs(rho1) + 1, c='black')
    ax.plot(rho1, np.ones_like(rho1) * -1, c='black')
    ax.plot(rho1, -(rho1**2 / 4), c='black')

    # Turn normal axes off
    for spine in ['left', 'bottom', 'top', 'right']:
        ax.spines[spine].set_visible(False)

    # Add arrows to represent axes
    axes_arrows = {'arrowstyle': '<|-|>', 'lw': 1.3}
    ax.annotate('', xy=(xmin, 0), xytext=(xmax, 0), arrowprops=axes_arrows)
    ax.annotate('', xy=(0, ymin), xytext=(0, ymax), arrowprops=axes_arrows)

    # Annotate the plot with equations
    plot_arrowsl = {'arrowstyle': '-|>', 'connectionstyle': "arc3, rad=-0.2"}
    plot_arrowsr = {'arrowstyle': '-|>', 'connectionstyle': "arc3, rad=0.2"}
    ax.annotate(r'$\rho_1 + \rho_2 < 1$', xy=(0.5, 0.3), xytext=(0.8, 0.6),
                arrowprops=plot_arrowsr, fontsize='12')
    ax.annotate(r'$\rho_1 + \rho_2 = 1$', xy=(0.38, 0.6), xytext=(0.6, 0.8),
                arrowprops=plot_arrowsr, fontsize='12')
    ax.annotate(r'$\rho_2 < 1 + \rho_1$', xy=(-0.5, 0.3), xytext=(-1.3, 0.6),
                arrowprops=plot_arrowsl, fontsize='12')
    ax.annotate(r'$\rho_2 = 1 + \rho_1$', xy=(-0.38, 0.6), xytext=(-1, 0.8),
                arrowprops=plot_arrowsl, fontsize='12')
    ax.annotate(r'$\rho_2 = -1$', xy=(1.5, -1), xytext=(1.8, -1.3),
                arrowprops=plot_arrowsl, fontsize='12')
    ax.annotate(r'${\rho_1}^2 + 4\rho_2 = 0$', xy=(1.15, -0.35),
                xytext=(1.5, -0.3), arrowprops=plot_arrowsr, fontsize='12')
    ax.annotate(r'${\rho_1}^2 + 4\rho_2 < 0$', xy=(1.4, -0.7),
                xytext=(1.8, -0.6), arrowprops=plot_arrowsr, fontsize='12')

    # Label categories of solutions
    ax.text(1.5, 1, 'Explosive\n growth', ha='center', fontsize=16)
    ax.text(-1.5, 1, 'Explosive\n oscillations', ha='center', fontsize=16)
    ax.text(0.05, -1.5, 'Explosive oscillations', ha='center', fontsize=16)
    ax.text(0.09, -0.5, 'Damped oscillations', ha='center', fontsize=16)

    # Add small marker to y-axis
    ax.axhline(y=1.005, xmin=0.495, xmax=0.505, c='black')
    ax.text(-0.12, -1.12, '-1', fontsize=10)
    ax.text(-0.12, 0.98, '1', fontsize=10)

    return fig

param_plot()
plt.show()
```

### Explanation of the graph

The graph portrays regions in which the $ (\lambda_1, \lambda_2) $
root pairs implied by the $ (\rho_1 = (a+b), \rho_2 = - b) $
difference equation parameter pairs in the Samuelson model are such that:

- $ (\lambda_1, \lambda_2) $ are complex with modulus less than
  $ 1 $ - in this case, the $ \{Y_t\} $ sequence displays damped
  oscillations
- $ (\lambda_1, \lambda_2) $ are both real, but one is strictly
  greater than $ 1 $ - this leads to explosive growth
- $ (\lambda_1, \lambda_2) $ are both real, but one is strictly
  less than $ -1 $ - this leads to explosive oscillations
- $ (\lambda_1, \lambda_2) $ are both real and both are less than
  $ 1 $ in absolute value - in this case, there is smooth
  convergence to the steady state without damped cycles

Later we’ll present the graph with a red mark showing the particular
point implied by the setting of $ (a,b) $

### Function to describe implications of characteristic polynomial

```python
def categorize_solution(rho1, rho2):
    """this function takes values of rho1 and rho2 and uses them to classify the type of solution"""

    discriminant = rho1 ** 2 + 4 * rho2
    if rho2 > 1 + rho1 or rho2 < -1:
        print('Explosive oscillations')
    elif rho1 + rho2 > 1:
        print('Explosive growth')
    elif discriminant < 0:
        print('Roots are complex with modulus less than one; therefore damped oscillations')
    else:
        print('Roots are real and absolute values are less than zero; therfore get smooth convergence to a steady state')
```

```python
### Test the categorize_solution function

categorize_solution(1.3, -.4)
```

    Roots are real and absolute values are less than zero; therfore get smooth convergence to a steady state

```none
Roots are real and absolute values are less than zero; therfore get smooth convergence to a steady state
```

### Function for plotting $ Y_t $ paths

A useful function for our work below

```python
def plot_y(function=None):
    """function plots path of Y_t"""
    plt.subplots(figsize=(12, 8))
    plt.plot(function)
    plt.xlabel('Time $t$')
    plt.ylabel('$Y_t$', rotation=0)
    plt.grid()
    plt.show()
```

### Manual or “by hand” root calculations

The following function calculates roots of the characteristic polynomial
using high school algebra

(We’ll calculate the roots in other ways later)

The function also plots a $ Y_t $ starting from initial conditions
that we set

```python
from cmath import sqrt

##=== This is a 'manual' method ===#

def y_nonstochastic(y_0=100, y_1=80, alpha=.92, beta=.5, gamma=10, n=80):

    """Takes values of parameters and computes roots of characteristic polynomial.
       It tells whether they are real or complex and whether they are less than unity in absolute value.
       It also computes a simulation of length n starting from the two given initial conditions for national income"""

    roots = []

    rho1 = alpha + beta
    rho2 = -beta

    print('rho_1 is ', rho1)
    print('rho_2 is ', rho2)

    discriminant = rho1 ** 2 + 4 * rho2

    if discriminant == 0:
        roots.append(-rho1 / 2)
        print('Single real root: ')
        print(''.join(str(roots)))
    elif discriminant > 0:
        roots.append((-rho1 + sqrt(discriminant).real) / 2)
        roots.append((-rho1 - sqrt(discriminant).real) / 2)
        print('Two real roots: ')
        print(''.join(str(roots)))
    else:
        roots.append((-rho1 + sqrt(discriminant)) / 2)
        roots.append((-rho1 - sqrt(discriminant)) / 2)
        print('Two complex roots: ')
        print(''.join(str(roots)))

    if all(abs(root) < 1 for root in roots):
        print('Absolute values of roots are less than one')
    else:
        print('Absolute values of roots are not less than one')

    def transition(x, t): return rho1 * x[t - 1] + rho2 * x[t - 2] + gamma

    y_t = [y_0, y_1]

    for t in range(2, n):
        y_t.append(transition(y_t, t))

    return y_t

plot_y(y_nonstochastic())
```

```none
rho_1 is  1.42
rho_2 is  -0.5
Two real roots:
[-0.6459687576256715, -0.7740312423743284]
Absolute values of roots are less than one
```

### Reverse engineering parameters to generate damped cycles

The next cell writes code that takes as inputs the modulus $ r $ and
phase $ \phi $ of a conjugate pair of complex numbers in polar form

$$
\lambda_1 = r \exp(i \phi), \quad \lambda_2 = r \exp(- i \phi)
$$

- The code assumes that these two complex numbers are the roots of the
  characteristic polynomial
- It then reverse engineers $ (a,b) $ and $ (\rho_1, \rho_2) $,
  pairs that would generate those roots

```python
### code to reverse engineer a  cycle
### y_t = r^t (c_1 cos(phi t) + c2 sin(phi t))
###

import cmath
import math

def f(r, phi):
    """
    Takes modulus r and angle phi of complex number r exp(j phi)
    and creates rho1 and rho2 of characteristic polynomial for which
    r exp(j phi) and r exp(- j phi) are complex roots.

    Returns the multiplier coefficient a and the accelerator coefficient b
    that verifies those roots.
    """
    g1 = cmath.rect(r, phi)  # Generate two complex roots
    g2 = cmath.rect(r, -phi)
    rho1 = g1 + g2           # Implied rho1, rho2
    rho2 = -g1 * g2
    b = -rho2                # Reverse engineer a and b that validate these
    a = rho1 - b
    return rho1, rho2, a, b

## Now let's use the function in an example
## Here are the example paramters

r = .95
period = 10   # Length of cycle in units of time
phi = 2 * math.pi/period

## Apply the function

rho1, rho2, a, b = f(r, phi)

print("a, b = ", a, b)
print("rho1, rho2 =", rho1, rho2)
```

    a, b =  (0.6346322893124001+0j) (0.9024999999999999-0j)
    rho1, rho2 = (1.5371322893124+0j) (-0.9024999999999999+0j)

```none
a, b =  (0.6346322893124001+0j) (0.9024999999999999-0j)
rho1, rho2 = (1.5371322893124+0j) (-0.9024999999999999+0j)
```

```python
## Print the real components of rho1 and rho2

rho1 = rho1.real
rho2 = rho2.real

rho1, rho2
```

```none
(1.5371322893124, -0.9024999999999999)
```

### Root finding using numpy

Here we’ll use numpy to compute the roots of the characteristic
polynomial

```python
r1, r2 = np.roots([1, -rho1, -rho2])

p1 = cmath.polar(r1)
p2 = cmath.polar(r2)

print("r, phi =", r, phi)
print("p1, p2 = ", p1, p2)
# print("g1, g2 = ", g1, g2)

print("a, b =", a, b)
print("rho1, rho2 =", rho1, rho2)
```

    r, phi = 0.95 0.6283185307179586
    p1, p2 =  (0.95, 0.6283185307179586) (0.95, -0.6283185307179586)
    a, b = (0.6346322893124001+0j) (0.9024999999999999-0j)
    rho1, rho2 = 1.5371322893124 -0.9024999999999999

```none
r, phi = 0.95 0.6283185307179586
p1, p2 =  (0.95, 0.6283185307179586) (0.95, -0.6283185307179586)
a, b = (0.6346322893124001+0j) (0.9024999999999999-0j)
rho1, rho2 = 1.5371322893124 -0.9024999999999999
```

```python
##=== This method uses numpy to calculate roots ===#

def y_nonstochastic(y_0=100, y_1=80, alpha=.9, beta=.8, gamma=10, n=80):

    """ Rather than computing the roots of the characteristic polynomial by hand as we did earlier, this function
    enlists numpy to do the work for us """

    # Useful constants
    rho1 = alpha + beta
    rho2 = -beta

    categorize_solution(rho1, rho2)

    # Find roots of polynomial
    roots = np.roots([1, -rho1, -rho2])
    print('Roots are', roots)

    # Check if real or complex
    if all(isinstance(root, complex) for root in roots):
        print('Roots are complex')
    else:
        print('Roots are real')

    # Check if roots are less than one
    if all(abs(root) < 1 for root in roots):
        print('Roots are less than one')
    else:
        print('Roots are not less than one')

    # Define transition equation
    def transition(x, t): return rho1 * x[t - 1] + rho2 * x[t - 2] + gamma

    # Set initial conditions
    y_t = [y_0, y_1]

    # Generate y_t series
    for t in range(2, n):
        y_t.append(transition(y_t, t))

    return y_t

plot_y(y_nonstochastic())
```

```none
Roots are complex with modulus less than one; therefore damped oscillations
Roots are [ 0.85+0.27838822j  0.85-0.27838822j]
Roots are complex
Roots are less than one
```

### Reverse engineered complex roots: example

The next cell studies the implications of revese engineered complex
roots

We’ll generate an **undamped** cycle of period 10

```python
r = 1   # generates undamped, nonexplosive cycles

period = 10   #  length of cycle in units of time
phi = 2 * math.pi/period

## Apply the reverse engineering function f

rho1, rho2, a, b = f(r, phi)

a = a.real  # drop the imaginary part so that it is a valid input into y_nonstochastic
b = b.real

print("a, b =", a, b)

ytemp = y_nonstochastic(alpha=a, beta=b, y_0=20, y_1=30)
plot_y(ytemp)
```

```none
a, b = 0.6180339887498949 1.0
Roots are complex with modulus less than one; therefore damped oscillations
Roots are [ 0.80901699+0.58778525j  0.80901699-0.58778525j]
Roots are complex
Roots are less than one
```

### Interactive demonstration

We'll use some widgets to show some of the consequences of parameters $(a,b)$ that are associated with complex conjugate roots of the characteristic polynomial. The widget will allow us to specify the amplitude and angle of the complex roots of the characteristic polynomial. (It will reverse engineer the $(a,b)$ that imply those roots.)

```python
from ipywidgets import interact

def choose_r_phi(r, phi):

    rho1, rho2, a, b = f(r, phi)
    a = a.real  # drop the imaginary part so that it is a valid input into y_nonstochastic
    b = b.real
    ytemp = y_nonstochastic(alpha=a, beta=b, y_0=20, y_1=30)

    return plot_y(ytemp)
```

```python
periodhighest = 20

interact(choose_r_phi, r=(0, 1.5), phi=(math.pi/periodhighest, 10 * math.pi/periodhighest))
```

    A Jupyter Widget

    <function __main__.choose_r_phi>

### Digression: using sympy to find roots

We can also use sympy to compute analytic formulas for the roots

```python
import sympy
from sympy import Symbol, init_printing
init_printing()

r1 = Symbol("rho_1")
r2 = Symbol("rho_2")
z = Symbol("z")

sympy.solve(z**2 - r1*z - r2, z)
```

$$
\left [ \frac{\rho_{1}}{2} - \frac{1}{2} \sqrt{\rho_{1}^{2} + 4 \rho_{2}},
\quad \frac{\rho_{1}}{2} + \frac{1}{2} \sqrt{\rho_{1}^{2} + 4 \rho_{2}}\right ]
$$

```python
a = Symbol("alpha")
b = Symbol("beta")
r1 = a + b
r2 = -b

sympy.solve(z**2 - r1*z - r2, z)
```

$$
\left [ \frac{\alpha}{2} + \frac{\beta}{2} - \frac{1}{2} \sqrt{\alpha^{2} +
2 \alpha \beta + \beta^{2} - 4 \beta}, \quad \frac{\alpha}{2} +
\frac{\beta}{2} + \frac{1}{2} \sqrt{\alpha^{2} + 2 \alpha \beta +
\beta^{2} - 4 \beta}\right ]
$$

## Stochastic shocks

Now we’ll construct some code to simulate the stochastic version of the
model that emerges when we add a random shock process to aggregate
demand

```python
def y_stochastic(y_0=0, y_1=0, alpha=0.8, beta=0.2, gamma=10, n=100, sigma=5):

    """This function takes parameters of a stochastic version of the model and proceeds to analyze
    the roots of the characteristic polynomial and also generate a simulation"""

    # Useful constants
    rho1 = alpha + beta
    rho2 = -beta

    # Categorize solution
    categorize_solution(rho1, rho2)

    # Find roots of polynomial
    roots = np.roots([1, -rho1, -rho2])
    print(roots)

    # Check if real or complex
    if all(isinstance(root, complex) for root in roots):
        print('Roots are complex')
    else:
        print('Roots are real')

    # Check if roots are less than one
    if all(abs(root) < 1 for root in roots):
        print('Roots are less than one')
    else:
        print('Roots are not less than one')

    # Generate shocks
    epsilon = np.random.normal(0, 1, n)

    # Define transition equation
    def transition(x, t): return rho1 * \
        x[t - 1] + rho2 * x[t - 2] + gamma + sigma * epsilon[t]

    # Set initial conditions
    y_t = [y_0, y_1]

    # Generate y_t series
    for t in range(2, n):
        y_t.append(transition(y_t, t))

    return y_t

plot_y(y_stochastic())
```

```none
Roots are real and absolute values are less than zero; therfore get smooth convergence to a steady state
[ 0.7236068  0.2763932]
Roots are real
Roots are less than one
```

Let’s do a simulation in which there are shocks and the characteristic polynomial has complex roots

```python
r = .97

period = 10   #  length of cycle in units of time
phi = 2 * math.pi/period

### apply the  reverse engineering function f

rho1, rho2, a, b = f(r, phi)

a = a.real  # drop the imaginary part so that it is a valid input into y_nonstochastic
b = b.real

print("a, b = ", a, b)
plot_y(y_stochastic(y_0=40, y_1 = 42, alpha=a, beta=b, sigma=2, n=100))
```

```none
a, b =  0.6285929690873979 0.9409000000000001
Roots are complex with modulus less than one; therefore damped oscillations
[ 0.78474648+0.57015169j  0.78474648-0.57015169j]
Roots are complex
Roots are less than one
```

## Government spending

This function computes a response to either a permanent or one-off increase
in government expenditures

```python
def y_stochastic_g(y_0=20,
                   y_1=20,
                   alpha=0.8,
                   beta=0.2,
                   gamma=10,
                   n=100,
                   sigma=2,
                   g=0,
                   g_t=0,
                   duration='permanent'):

    """This program computes a response to a permanent increase in government expenditures that occurs
       at time 20"""

    # Useful constants
    rho1 = alpha + beta
    rho2 = -beta

    # Categorize solution
    categorize_solution(rho1, rho2)

    # Find roots of polynomial
    roots = np.roots([1, -rho1, -rho2])
    print(roots)

    # Check if real or complex
    if all(isinstance(root, complex) for root in roots):
        print('Roots are complex')
    else:
        print('Roots are real')

    # Check if roots are less than one
    if all(abs(root) < 1 for root in roots):
        print('Roots are less than one')
    else:
        print('Roots are not less than one')

    # Generate shocks
    epsilon = np.random.normal(0, 1, n)

    def transition(x, t, g):

        # Non-stochastic - separated to avoid generating random series when not needed
        if sigma == 0:
            return rho1 * x[t - 1] + rho2 * x[t - 2] + gamma + g

        # Stochastic
        else:
            epsilon = np.random.normal(0, 1, n)
            return rho1 * x[t - 1] + rho2 * x[t - 2] + gamma + g + sigma * epsilon[t]

    # Create list and set initial conditions
    y_t = [y_0, y_1]

    # Generate y_t series
    for t in range(2, n):

        # No government spending
        if g == 0:
            y_t.append(transition(y_t, t))

        # Government spending (no shock)
        elif g != 0 and duration == None:
            y_t.append(transition(y_t, t))

        # Permanent government spending shock
        elif duration == 'permanent':
            if t < g_t:
                y_t.append(transition(y_t, t, g=0))
            else:
                y_t.append(transition(y_t, t, g=g))

        # One-off government spending shock
        elif duration == 'one-off':
            if t == g_t:
                y_t.append(transition(y_t, t, g=g))
            else:
                y_t.append(transition(y_t, t, g=0))
    return y_t
```

A permanent government spending shock can be simulated as follows

```python
plot_y(y_stochastic_g(g=10, g_t=20, duration='permanent'))
```

```none
Roots are real and absolute values are less than zero; therfore get smooth convergence to a steady state
[ 0.7236068  0.2763932]
Roots are real
Roots are less than one
```

We can also see the response to a one time jump in government expenditures

```python
plot_y(y_stochastic_g(g=500, g_t=50, duration='one-off'))
```

```none
Roots are real and absolute values are less than zero; therfore get smooth convergence to a steady state
[ 0.7236068  0.2763932]
Roots are real
Roots are less than one
```

## Wrapping everything into a class

Up to now we have written functions to do the work

Now we’ll roll up our sleeves and write a Python class called Samuelson
for the Samuleson model

```python
class Samuelson():

    r"""This class represents the Samuelson model, otherwise known as the
    multiple-accelerator model. The model combines the Keynesian multipler
    with the accelerator theory of investment.

    The path of output is governed by a linear second-order difference equation

    .. math::

        Y_t =  + \alpha (1 + \beta) Y_{t-1} - \alpha \beta Y_{t-2}

    Parameters
    ----------
    y_0 : scalar
        Initial condition for Y_0
    y_1 : scalar
        Initial condition for Y_1
    alpha : scalar
        Marginal propensity to consume
    beta : scalar
        Accelerator coefficient
    n : int
        Number of iterations
    sigma : scalar
        Volatility parameter. Must be greater than or equal to 0. Set
        equal to 0 for non-stochastic model.
    g : scalar
        Government spending shock
    g_t : int
        Time at which government spending shock occurs. Must be specified
        when duration != None.
    duration : {None, 'permanent', 'one-off'}
        Specifies type of government spending shock. If none, government
        spending equal to g for all t.

    """

    def __init__(self,
                 y_0=100,
                 y_1=50,
                 alpha=1.3,
                 beta=0.2,
                 gamma=10,
                 n=100,
                 sigma=0,
                 g=0,
                 g_t=0,
                 duration=None):

        self.y_0, self.y_1, self.alpha, self.beta = y_0, y_1, alpha, beta
        self.n, self.g, self.g_t, self.duration = n, g, g_t, duration
        self.gamma, self.sigma = gamma, sigma
        self.rho1 = alpha + beta
        self.rho2 = -beta
        self.roots = np.roots([1, -self.rho1, -self.rho2])

    def root_type(self):
        if all(isinstance(root, complex) for root in self.roots):
            return 'Complex conjugate'
        elif len(self.roots) > 1:
            return 'Double real'
        else:
            return 'Single real'

    def root_less_than_one(self):
        if all(abs(root) < 1 for root in self.roots):
            return True

    def solution_type(self):
        rho1, rho2 = self.rho1, self.rho2
        discriminant = rho1 ** 2 + 4 * rho2
        if rho2 >= 1 + rho1 or rho2 <= -1:
            return 'Explosive oscillations'
        elif rho1 + rho2 >= 1:
            return 'Explosive growth'
        elif discriminant < 0:
            return 'Damped oscillations'
        else:
            return 'Steady state'

    def _transition(self, x, t, g):

        # Non-stochastic - separated to avoid generating random series when not needed
        if self.sigma == 0:
            return self.rho1 * x[t - 1] + self.rho2 * x[t - 2] + self.gamma + g

        # Stochastic
        else:
            epsilon = np.random.normal(0, 1, self.n)
            return self.rho1 * x[t - 1] + self.rho2 * x[t - 2] + self.gamma + g + self.sigma * epsilon[t]

    def generate_series(self):

        # Create list and set initial conditions
        y_t = [self.y_0, self.y_1]

        # Generate y_t series
        for t in range(2, self.n):

            # No government spending
            if self.g == 0:
                y_t.append(self._transition(y_t, t))

            # Government spending (no shock)
            elif self.g != 0 and self.duration == None:
                y_t.append(self._transition(y_t, t))

            # Permanent government spending shock
            elif self.duration == 'permanent':
                if t < self.g_t:
                    y_t.append(self._transition(y_t, t, g=0))
                else:
                    y_t.append(self._transition(y_t, t, g=self.g))

            # One-off government spending shock
            elif self.duration == 'one-off':
                if t == self.g_t:
                    y_t.append(self._transition(y_t, t, g=self.g))
                else:
                    y_t.append(self._transition(y_t, t, g=0))
        return y_t

    def summary(self):
        print('Summary\n' + '-' * 50)
        print('Root type: ' + self.root_type())
        print('Solution type: ' + self.solution_type())
        print('Roots: ' + str(self.roots))

        if self.root_less_than_one() == True:
            print('Absolute value of roots is less than one')
        else:
            print('Absolute value of roots is not less than one')

        if self.sigma > 0:
            print('Stochastic series with sigma = ' + str(self.sigma))
        else:
            print('Non-stochastic series')

        if self.g != 0:
            print('Government spending equal to ' + str(self.g))

        if self.duration != None:
            print(self.duration.capitalize() +
                  ' government spending shock at t = ' + str(self.g_t))

    def plot(self):
        fig, ax = plt.subplots(figsize=(12, 8))
        ax.plot(self.generate_series())
        ax.set(xlabel='Iteration', xlim=(0, self.n))
        ax.set_ylabel('$Y_t$', rotation=0)
        ax.grid()

        # Add parameter values to plot
        paramstr = '$\\alpha=%.2f$\n$\\beta=%.2f$\n$\\gamma=%.2f$\n$\\sigma=%.2f$\n$\\rho_1=%.2f$\n$\\rho_2=%.2f$'%(
                    self.alpha, self.beta, self.gamma, self.sigma, self.rho1, self.rho2)
        props = dict(fc='white', pad=10, alpha=0.5)
        ax.text(0.87, 0.05, paramstr, transform=ax.transAxes,
                fontsize=12, bbox=props, va='bottom')

        return fig

    def param_plot(self):

        # Uses the param_plot() function defined earlier (it is then able
        # to be used standalone or as part of the model)

        fig = param_plot()
        ax = fig.gca()

        # Add lambda values to legend
        for i, root in enumerate(self.roots):
            if isinstance(root, complex):
                operator = ['+', ''] # Need to fill operator for positive as string is split apart
                label = r'$\lambda_{0} = {1.real:.2f} {2} {1.imag:.2f}i$'.format(i+1, sam.roots[i], operator[i])
            else:
                label = r'$\lambda_{0} = {1.real:.2f}$'.format(i+1, sam.roots[i])
            ax.scatter(0, 0, 0, label=label) # dummy to add to legend

        # Add rho pair to plot
        ax.scatter(self.rho1, self.rho2, 100, 'red', '+', label=r'$(\ \rho_1, \ \rho_2 \ )$', zorder=5)

        plt.legend(fontsize=12, loc=3)

        return fig
```

### Illustration of Samuelson class

Now we’ll put our Samuelson class to work on an example

```python
sam = Samuelson(alpha=0.8, beta=0.5, sigma=2, g=10, g_t=20, duration='permanent')
sam.summary()
```

    Summary
    --------------------------------------------------
    Root type: Complex conjugate
    Solution type: Damped oscillations
    Roots: [ 0.65+0.27838822j  0.65-0.27838822j]
    Absolute value of roots is less than one
    Stochastic series with sigma = 2
    Government spending equal to 10
    Permanent government spending shock at t = 20

```none
Summary
---------------------------------------------------------
Root type: Complex conjugate
Solution type: Damped oscillations
Roots: [ 0.65+0.27838822j  0.65-0.27838822j]
Absolute value of roots is less than one
Stochastic series with sigma = 2
Government spending equal to 10
Permanent government spending shock at t = 20
```

```python
sam.plot()
plt.show()
```

### Using the graph

We’ll use our graph to show where the roots lie and how their location
is consistent with the behavior of the path just graphed

The red $+$ sign shows the location of the roots

```python
sam.param_plot()
plt.show()
```

## Using the LinearStateSpace class

It turns out that we can use the [QuantEcon.py](http://quantecon.org/python_index.html)
[LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class to do
much of the work that we have done from scratch above

Here is how we map the Samuelson model into an instance of a
LinearStateSpace class

```python
from quantecon import LinearStateSpace

""" This script maps the Samuelson model in the the LinearStateSpace class"""
alpha = 0.8
beta = 0.9
rho1 = alpha + beta
rho2 = -beta
gamma = 10
sigma = 1
g = 10
n = 100

A = [[1,        0,        0],
     [gamma + g, rho1, rho2],
     [0,        1,        0]]

G = [[gamma + g, rho1, rho2],         # this is Y_{t+1}
     [gamma,    alpha,    0],         # this is C_{t+1}
     [0,     beta,    -beta]]         # this is I_{t+1}

mu_0 = [1, 100, 100]
C = np.zeros((3,1))
C[1] = sigma # stochastic

sam_t = LinearStateSpace(A, C, G, mu_0=mu_0)

x, y = sam_t.simulate(ts_length=n)

fig, axes = plt.subplots(3, 1, sharex=True, figsize=(15, 8))
titles = ['Output ($Y_t$)', 'Consumption ($C_t$)', 'Investment ($I_t$)']
colors = ['darkblue', 'red', 'purple']
for ax, series, title, color in zip(axes, y, titles, colors):
    ax.plot(series, color=color)
    ax.set(title=title, xlim=(0, n))
    ax.grid()

axes[-1].set_xlabel('Iteration')

plt.show()
```

### Other methods in the LinearStateSpace class

Let’s plot **impulse response functions** for the instance of the
Samuelson model using a method in the LinearStateSpace class

```python
imres = sam_t.impulse_response()
imres = np.asarray(imres)
y1 = imres[:, :, 0]
y2 = imres[:, :, 1]
y1.shape
```

$$
\left ( 2, \quad 6, \quad 1\right )
$$

Now let’s compute the zeros of the characteristic polynomial by simply
calculating the eigenvalues of $ A $

```python
A = np.asarray(A)
w, v = np.linalg.eig(A)
print(w)
```

    [ 0.85+0.42130749j  0.85-0.42130749j  1.00+0.j        ]

```none
[ 0.85+0.42130749j  0.85-0.42130749j  1.00+0.j        ]
```

### Inheriting methods from LinearStateSpace

We could also create a subclass of LinearStateSpace (inheriting all its
methods and attributes) to add more functions to use

```python
class SamuelsonLSS(LinearStateSpace):

    """
    this subclass creates a Samuelson multiplier-accelerator model
    as a linear state space system
    """
    def __init__(self,
                 y_0=100,
                 y_1=100,
                 alpha=0.8,
                 beta=0.9,
                 gamma=10,
                 sigma=1,
                 g=10):

        self.alpha, self.beta = alpha, beta
        self.y_0, self.y_1, self.g = y_0, y_1, g
        self.gamma, self.sigma = gamma, sigma

        # Define intial conditions
        self.mu_0 = [1, y_0, y_1]

        self.rho1 = alpha + beta
        self.rho2 = -beta

        # Define transition matrix
        self.A = [[1,                 0,         0],
                  [gamma + g, self.rho1, self.rho2],
                  [0,                 1,         0]]

        # Define output matrix
        self.G = [[gamma + g, self.rho1, self.rho2],         # this is Y_{t+1}
                  [gamma,         alpha,         0],         # this is C_{t+1}
                  [0,              beta,     -beta]]         # this is I_{t+1}

        self.C = np.zeros((3, 1))
        self.C[1] = sigma  # stochastic

        # Initialize LSS with parameters from Samuleson model
        LinearStateSpace.__init__(self, self.A, self.C, self.G, mu_0=self.mu_0)

    def plot_simulation(self, ts_length=100, stationary=True):

        # Temporarily store original parameters
        temp_mu = self.mu_0
        temp_Sigma = self.Sigma_0

        # Set distribution parameters equal to their stationary values for simulation
        if stationary == True:
            try:
                self.mu_x, self.mu_y, self.sigma_x, self.sigma_y = self.stationary_distributions()
                self.mu_0 = self.mu_y
                self.Sigma_0 = self.sigma_y
            # Exception where no convergence achieved when calculating stationary distributions
            except ValueError:
                print('Stationary distribution does not exist')

        x, y = self.simulate(ts_length)

        fig, axes = plt.subplots(3, 1, sharex=True, figsize=(15, 8))
        titles = ['Output ($Y_t$)', 'Consumption ($C_t$)', 'Investment ($I_t$)']
        colors = ['darkblue', 'red', 'purple']
        for ax, series, title, color in zip(axes, y, titles, colors):
            ax.plot(series, color=color)
            ax.set(title=title, xlim=(0, n))
            ax.grid()

        axes[-1].set_xlabel('Iteration')

        # Reset distribution parameters to their initial values
        self.mu_0 = temp_mu
        self.Sigma_0 = temp_Sigma

        return fig

    def plot_irf(self, j=5):

        x, y = self.impulse_response(j)

        # Reshape into 3 x j matrix for plotting purposes
        yimf = np.array(y).flatten().reshape(j+1, 3).T

        fig, axes = plt.subplots(3, 1, sharex=True, figsize=(15, 8))
        labels = ['$Y_t$', '$C_t$', '$I_t$']
        colors = ['darkblue', 'red', 'purple']
        for ax, series, label, color in zip(axes, yimf, labels, colors):
            ax.plot(series, color=color)
            ax.set(xlim=(0, j))
            ax.set_ylabel(label, rotation=0, fontsize=14, labelpad=10)
            ax.grid()

        axes[0].set_title('Impulse Response Functions')
        axes[-1].set_xlabel('Iteration')

        return fig

    def multipliers(self, j=5):
        x, y = self.impulse_response(j)
        return np.sum(np.array(y).flatten().reshape(j+1, 3), axis=0)
```

### Illustrations

Let’s show how we can use the SamuelsonLSS

```python
samlss = SamuelsonLSS()
```

```python
samlss.plot_simulation(100, stationary=False)
plt.show()
```

```python
samlss.plot_simulation(100, stationary=True)
plt.show()
```

```python
samlss.plot_irf(100)
plt.show()
```

```python
samlss.multipliers()
```

    array([ 7.414389,  6.835896,  0.578493])

```none
array([ 7.414389,  6.835896,  0.578493])
```

## Pure multiplier model

Let’s shut down the accelerator by setting $ b=0 $ to get a pure
multiplier model

- the absence of cycles gives an idea about why Samuelson included the
  accelerator

```python
pure_multiplier = SamuelsonLSS(alpha=0.95, beta=0)
```

```python
pure_multiplier.plot_simulation()
```

```none
Stationary distribution does not exist
```

```python
pure_multiplier = SamuelsonLSS(alpha=0.8, beta=0)
```

```python
pure_multiplier.plot_simulation()
```

```python
pure_multiplier.plot_irf(100)
```

## Summary

In this lecture, we wrote functions and classes to represent non-stochastic and
stochastic versions of the Samuelson (1939) multiplier-accelerator model, described
in [[Sam39]](zreferences.ipynb#samuelson1939)

We saw that different parameter values led to different output paths, which
could either be stationary, explosive, or ocillating

We also were able to represent the model using the [QuantEcon.py](http://quantecon.org/python_index.html)
[LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class


text:
# Point particle on rotating curve

## Point particle on arbitrary curve

Let us assume that point with mass $m=1$ moves freely on a flat curve which spins around vertical axis with angular velocity $\omega_0$.
In general case we have a curve given by an implicit equation in a form:
$$ f(x,y) = 0 $$

We can treat the problem in 2d in the rotating frame of reference. However, we need to take into account additional force in non-inertial system which acts in $x$ direction:

$$F_x = m x \omega_0^2.$$

```python
load('cas_utils.sage')
```

```python
var('t')
var('g')
var('w0', latex_name='\omega_0')
xy_names = [('x','x'),('y','y')]
```

```python
to_fun, to_var = make_symbols(xy_names)
```

dAlembert principle takes the following form:

```python
dAlemb = (X.diff(t,2)-x*w0^2)*dx + (Y.diff(t,2)+g)*dy
```

```python
showmath(dAlemb.subs(to_var))
```

We will need differential form of contraints, which can also be obtained by computer algebra.
We will use implicit equation of a circle of unit radius. In this case we can compare solutions obtained in generic procedure with one computer in generalized coordinates.

```python
f = x^2+ (y-1)^2 - 1

plt_constraints = implicit_plot(f,(x,-1,1),(y,0,2),figsize=5,color='gray')
plt_constraints
```

We will use constaints to connect variations $\delta x$ and $\delta y$ with each other, therefore we use formula for complete differential:

```python
df = diff(f,x)*dx+diff(f,y)*dy
showmath(df)
```

We solve above equation for $\delta x$ and substitute it to dAlembert equation. Since there is only one variation left, then its coeficcient must be zero. It is a first differential equation.

```python
eq1 = dAlemb.subs(df.solve(dx)).coefficient(dy).subs(to_var)
showmath(eq1)
```

The second equation can be obtained, for example by differentiation of contraints in time two times.

```python
eq2 = f.subs(to_fun).diff(t,2).subs(to_var)
showmath(eq2)
```

We solve above system for $\ddot x$ and $\ddot y$:

```python
sol = solve( [eq1,eq2], [xdd,ydd])[0]
showmath( sol[0] )
```

```python
showmath( sol[1] )
```

We are ready to set numerical values of parameters and solve system of ODEs:

```python
times = srange(0,24,0.01)

ode = [xd,yd] + [ sol_.rhs().subs({w0:4.2,g:9.81}) for sol_ in sol]

x0 = 0.6
y0 = f.subs({x:x0}).solve(y)[1].rhs().n()
dx0 = 0.0
dy0 = df.subs({x:x0,y:y0,dx:dx0}).solve(dy)[0].rhs().n()

numsol = desolve_odeint(ode,[x0,y0,dx0,dy0], times, [x,y,xd,yd])
p = line(zip(numsol[:,0],numsol[:,1]), color='red')
p2 = line(zip(times,numsol[:,0]), figsize=5, ymin=-1, ymax=1)

(p+plt_constraints).show(title='phase space')
p2.show(title=r'$x(t)$')
```

**Experiment with Sage!**

Change parameters to see if other types of solution can be found!

### Explicit form of contraints

If contraits are in explicit for we can simplify the procedure.

```python
# f = 1/2*(x)^2 - y
f = 1-sqrt(1-x^2) - y
```

```python
df = diff(f,x)*dx+diff(f,y)*dy
```

```python
eq1 = dAlemb.subs(df.solve(dx)).coefficient(dy).subs(to_var)
showmath(eq1)
```

```python
eq2 = f.subs(to_fun).diff(t,2).subs(to_var)
showmath(eq2)
```

```python
sol = solve( [eq1,eq2], [xdd,ydd])[0]
showmath( sol[0] )
```

We see that the first equation does not depent on $y$. We can solve it numericaly.

```python
y_on_curve = y.subs(f.solve(y)[0])
```

```python
y_fun = fast_callable(y_on_curve, vars = [x])
```

```python
y_fun(numsol[:,0])
```

```python
ode = [xd,sol[0].rhs().subs({w0:4.2,g:9.81})]

times = srange(0,24,0.01)
numsol = desolve_odeint(ode,[.6,.0],times,[x,xd])
p = line(zip(numsol[:,0],y_fun(numsol[:,0])),color='red')
p2 = line(zip(times,numsol[:,0]),figsize=5,ymin=-1,ymax=1)

(p+plt_constraints).show(title='phase space')
p2.show(title=r'$x(t)$')

```

We see that it was possible to obtain the same solution. However, the implicit form allows to solve the system for  $y>1$. The explicit formulation will fail in wuch a case.

## Point particle on rotating circle

If particle moves on the rotating circle, we can use the symmetry and transform dAlembert principle to  polar coordinates in rotating frame.

```python
var('x y t')
var('w0 l g')

xy_names = ['x','y']
uv_names = [('phi','\phi')]

to_fun, to_var = make_symbols(xy_names,uv_names)
```

```python
x2u = {x:l*sin(phi),y:1-l*cos(phi)}
```

```python
showmath(x2u)
```

```python
transform_virtual_displacements(xy_names,uv_names,verbose=True)
```

```python
dAlemb = (x.subs(x2u).subs(to_fun).diff(t,2)-w0^2*x.subs(x2u))*dx_polar + \
         (y.subs(x2u).subs(to_fun).diff(t,2)+g)*dy_polar
dAlemb = dAlemb.subs(to_var)
show(dAlemb)
```

```python
sol = dAlemb.expand().coefficient(dphi).trig_simplify().solve(phidd)
show( sol)
```

### Efffective potential

In polar coordinates we can analyze the system without solving. Namely we can observe that the generalized force depends only on $\phi$, thus:
\begin{equation}
{\ddot \phi} = \underbrace{\frac{l {\omega_0}^{2} \cos\left({\phi}\right) \sin\left({\phi}\right) - g \sin\left({\phi}\right)}{l}}_{-\displaystyle\frac{\partial U_{eff}(\phi)}{\partial\phi}}
\end{equation}

we can interpret the equation of motion as a motion in an effective potential.

Let's extrat this potential and see how it behaves for different parameters:

```python
Ueff = -sol[0].rhs().expand().integrate(phi)
showmath(Ueff)
```

```python
Ueff = -sol[0].rhs().expand().subs({g:9.81,l:1}).integrate(phi)
showmath(Ueff)
```

```python
plot( [Ueff.subs(w0==w0_) for w0_ in [0,4,8,12]], (phi,-pi,pi),\
     legend_label=[0,4,8,12], figsize=(6,3),\
     tick_formatter=[pi,None],ticks=[pi/2,None])
```

### Numerical solutions

```python
arcsin(0.4)
```

```python
times = srange(0,24,0.01)
pars = {l:1, w0:4.2, g:9.81}

x_phi = fast_callable( x.subs(x2u).subs(pars),vars=[phi])
y_phi = fast_callable( y.subs(x2u).subs(pars),vars=[phi])

ode = [phid,sol[0].rhs().subs(pars)]
show(ode)
numsol = desolve_odeint(ode,[arcsin(0.6),0],times,[phi,phid])
p = line(zip(x_phi(numsol[:,0]),y_phi(numsol[:,0])),color='red')
p2 = line(zip(times,x_phi(numsol[:,0])),figsize=4,ymin=-pi/2,ymax=pi/2)

(p+plt_constraints).show()
p2.show()
```

## Lagrange approach

### Rotating system of coordinates

In a rotating system of coordinates, the bead moves on the circle in $(x,y)$ plane. At the same time there is a rotation with angular velocity $\omega_0$ around vetical axis. It means that a bead will have addtional velocity, perpedicular to $(x,y)$ plane, which magnitute is:
$$
v_{rot} = x\omega_0 = l {\omega_0} \sin\left({\phi}\right)
$$

```python
Ek = 1/2*sum([x_.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 for x_ in [x,y]])
Ek += 1/2*(x.subs(x2u)*w0)^2
Ek = Ek.trig_simplify()
showmath(Ek)
```

```python
Ep = g*y.subs(x2u)
showmath(Ep)
```

```python
L = Ek-Ep
showmath(L)
```

```python
EL = L.diff(phid).subs(to_fun).diff(t).subs(to_var) - L.diff(phi)
```

```python
showmath(EL)
```

### Code generation

We can readily generate code which can be used in external programs.

En example can be found in [simulation and 3d vis](beadk3d.ipynb)

```python
from sympy import ccode
oderhs = sol[0].rhs()
ccode(oderhs._sympy_())
```

\newpage


text:
# Solution of the problem 6-15 of Michael Rider book, using the methods of dynamics of rigid bodies.
# by **Oscar Alberto Santos Muñoz**

This code is licensed by MIT LICENSE

Copyright (c) 2020 OscarSantos98

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

    Run in Google Colab

    View source on GitHub

```
import sympy as sp
```

From the problem we know that $\omega2 = 9 \; \frac{rad}{s} \; c.c.w.$ and therefore $\alpha2 = 0 \; \frac{rad}{s^{2}} $

# In the next line the values of the magnitude of the vectors are introduced.

```
a2, a3, a4, a5, a6, b3 = 70, 100, 80, 120, 100, 65
```

Print the values of the magnitudes.

```
print("Links length of the mechanism\na2 =",a2,"\na3 =",a3,"\na4 =",a4,"\na5 =",a5,"\na6 =",a6,"\nb3 =",b3)
```

    Links length of the mechanism
    a2 = 70
    a3 = 100
    a4 = 80
    a5 = 120
    a6 = 100
    b3 = 65

```
theta2 = (95*sp.pi/180).evalf()
theta3 = (4.63498*sp.pi/180).evalf()
theta4 = (76.57626*sp.pi/180).evalf()
theta5 = (-7.2327*sp.pi/180).evalf()
theta6 = (105.74221*sp.pi/180).evalf()
beta = (35.18382*sp.pi/180).evalf()
```

Create alpha, beta and theta strings.

```
th_str = []
w_str = []
alp_str = []
for i in range(2,7):
	th_str.append(chr(952) + str(i))
	w_str.append(chr(969) + str(i))
	alp_str.append(chr(945) + str(i))
```

Print the values of the angles.

```
print("Angles\n",th_str[0]," =",theta2,"\n",th_str[1]," =",theta3,"\n",th_str[2]," =",theta4,"\n",th_str[3]," =",theta5,"\n",th_str[4]," =",theta6,"\n",chr(946)," =",beta)
```

    Angles
     θ2  = 1.65806278939461
     θ3  = 0.0808956617640868
     θ4  = 1.33650786586323
     θ5  = -0.126234428808994
     θ6  = 1.84554972283527
     β  = 0.614073502429031

# The necessary position vectors will be computed.

```
rCO21 = sp.Matrix(1,3,[a2*sp.cos(theta2),a2*sp.sin(theta2),0])
rDC = sp.Matrix(1,3,[a3*sp.cos(theta3),a3*sp.sin(theta3),0])
rDO41 = sp.Matrix(1,3,[a4*sp.cos(theta4),a4*sp.sin(theta4),0])
rPC = sp.Matrix(1,3,[b3*sp.cos(theta3+beta),b3*sp.sin(theta3+beta),0])
rFP = sp.Matrix(1,3,[a5*sp.cos(theta5),a5*sp.sin(theta5),0])
rFO61 = sp.Matrix(1,3,[a6*sp.cos(theta6),a6*sp.sin(theta6),0])
```

Print the results of the position vectors.

```
print("Position vectors\nrCO21 =",rCO21,"\nrDC =",rDC,"\nrDO41 =",rDO41,"\nrPC =",rPC,"\nrFP =",rFP,"\nrFO61 =",rFO61)
```

    Position vectors
    rCO21 = Matrix([[-6.10090199233608, 69.7336288664222, 0]])
    rDC = Matrix([[99.6729729954293, 8.08074589703349, 0]])
    rDO41 = Matrix([[18.5720755743718, 77.8143817610846, 0]])
    rPC = Matrix([[49.9247741348628, 41.6235141186197, 0]])
    rFP = Matrix([[119.045161106793, -15.1079322561889, 0]])
    rFO61 = Matrix([[-27.1309590837080, 96.2492132913208, 0]])

# Velocity analysis: Definition of the angular velocities.

Declare symbolic variables.

```
w2, w3, w4, w5, w6 = sp.symbols('w2 w3 w4 w5 w6')
```

Create vectors with the known values and the symbolic variables.

```
omega2 = sp.Matrix(1,3,[0,0,9])
omega3 = sp.Matrix(1,3,[0,0,w3])
omega4 = sp.Matrix(1,3,[0,0,w4])
omega5 = sp.Matrix(1,3,[0,0,w5])
omega6 = sp.Matrix(1,3,[0,0,w6])
```

Print the vectors.

```
print("Angular velocities\n",w_str[0],"=",omega2,"\n",w_str[1],"=",omega3,"\n",w_str[2],"=",omega4,"\n",w_str[3],"=",omega5,"\n",w_str[4],"=",omega6)
```

    Angular velocities
     ω2 = Matrix([[0, 0, 9]])
     ω3 = Matrix([[0, 0, w3]])
     ω4 = Matrix([[0, 0, w4]])
     ω5 = Matrix([[0, 0, w5]])
     ω6 = Matrix([[0, 0, w6]])

Calculate velocities.

```
VC = omega2.cross(rCO21)
VD1 = VC + omega3.cross(rDC)
VD2 = omega4.cross(rDO41)
```

Print the results of the point velocities

```
print("Point velocities\nVC = ",w_str[0]," x rCO21 = ",VC,"\nVD1 = VC + ",w_str[1]," x rDC = ",VD1,"\nVD2 = ",w_str[2]," x rDO41 = ",VD2)
```

    Point velocities
    VC =  ω2  x rCO21 =  Matrix([[-627.602659797800, -54.9081179310247, 0]])
    VD1 = VC +  ω3  x rDC =  Matrix([[-8.08074589703349*w3 - 627.6026597978, 99.6729729954293*w3 - 54.9081179310247, 0]])
    VD2 =  ω4  x rDO41 =  Matrix([[-77.8143817610846*w4, 18.5720755743718*w4, 0]])

Solve for w3 and w4.

```
Eqv1 = sp.Eq(VD1[0],VD2[0])
Eqv2 = sp.Eq(VD1[1],VD2[1])
answ = sp.solve([Eqv1,Eqv2])
```

Print the answer.

```
print("Solution for w3 and w4 of matching velocity equations of point D\n",answ)
```

    Solution for w3 and w4 of matching velocity equations of point D
     {w3: 2.09422879753408, w4: 8.28285949170457}

Substitue the answer in the vector w3 and w4.

```
omega3 = omega3.subs(w3,answ[w3])
omega4 = omega4.subs(w4,answ[w4])
```

```
print("Substituting w3 and w4 in vectors",w_str[1]," and ",w_str[2],"\n",w_str[1],"=",omega3,"\n",w_str[2],"=",omega4)
```

    Substituting w3 and w4 in vectors ω3  and  ω4
     ω3 = Matrix([[0, 0, 2.09422879753408]])
     ω4 = Matrix([[0, 0, 8.28285949170457]])

Calculate the rest of the velocities.

```
VP = VC + omega3.cross(rPC)
VF1 = VP + omega5.cross(rFP)
VF2 = omega6.cross(rFO61)
```

Print the result of the velocities.

```
print("Point velocities\nVP = VC + ",w_str[1]," x rPC = ",VP,"\nVF1 = VP + ",w_str[3]," x rFP = ",VF1,"\nVF2 = ",w_str[4]," x rFO61 = ",VF2)
```

    Point velocities
    VP = VC +  ω3  x rPC =  Matrix([[-714.771821719580, 49.6457817725896, 0]])
    VF1 = VP +  ω5  x rFP =  Matrix([[15.1079322561889*w5 - 714.77182171958, 119.045161106793*w5 + 49.6457817725896, 0]])
    VF2 =  ω6  x rFO61 =  Matrix([[-96.2492132913208*w6, -27.130959083708*w6, 0]])

Solve for w5 and w6.

```
Eqv3 = sp.Eq(VF1[0],VF2[0])
Eqv4 = sp.Eq(VF1[1],VF2[1])
answ2 = sp.solve([Eqv3,Eqv4])
```

Print the answer.

```
print("Solution for w5 and w6 of matching velocity equations of point F\n",answ2)
```

    Solution for w5 and w6 of matching velocity equations of point F
     {w5: -2.18777796287317, w6: 7.76966997860839}

Substitue the answer in the vector w5 and w5.

```
omega5 = omega5.subs(w5,answ2[w5])
omega6 = omega6.subs(w6,answ2[w6])
```

```
print("Substituting w5 and w6 in vectors",w_str[3]," and ",w_str[4],"\n",w_str[3],"=",omega5,"\n",w_str[4],"=",omega6)
```

    Substituting w5 and w6 in vectors ω5  and  ω6
     ω5 = Matrix([[0, 0, -2.18777796287317]])
     ω6 = Matrix([[0, 0, 7.76966997860839]])

# Acceleration analysis. Definition of the angular accelerations.

Declare symbolic variables.

```
aa2, aa3, aa4, aa5, aa6 = sp.symbols('aa2 aa3 aa4 aa5 aa6')
```

Create vectors with the known values and the symbolic variables.

```
alpha2 = sp.Matrix(1,3,[0,0,0])
alpha3 = sp.Matrix(1,3,[0,0,aa3])
alpha4 = sp.Matrix(1,3,[0,0,aa4])
alpha5 = sp.Matrix(1,3,[0,0,aa5])
alpha6 = sp.Matrix(1,3,[0,0,aa6])
```

Print the vectors.

```
print("Angular accelerations\n",alp_str[0],"=",alpha2,"\n",alp_str[1],"=",alpha3,"\n",alp_str[2],"=",alpha4,"\n",alp_str[3],"=",alpha5,"\n",alp_str[4],"=",alpha6)
```

    Angular accelerations
     α2 = Matrix([[0, 0, 0]])
     α3 = Matrix([[0, 0, aa3]])
     α4 = Matrix([[0, 0, aa4]])
     α5 = Matrix([[0, 0, aa5]])
     α6 = Matrix([[0, 0, aa6]])

Calculate accelerations.

```
AC = alpha2.cross(rCO21) + omega2.cross(omega2.cross(rCO21))
AD3 = AC + alpha3.cross(rDC) + omega3.cross(omega3.cross(rDC))
AD4 = alpha4.cross(rDO41) + omega4.cross(omega4.cross(rDO41))
```

Print the results of the accelerations.

```
print("Point accelerations\nAC = ",alp_str[0]," x rCO21 + ",w_str[0]," x ",w_str[0]," x rCO21 = ",AC,"\nAD3 = AC + ",alp_str[1]," x rDC + ",w_str[1]," x ",w_str[1]," x rDC = ",AD3,"\nAD4 = ",alp_str[2]," x rDO41 + ",w_str[2]," x ",w_str[2]," x rDO41 = ",AD4)
```

    Point accelerations
    AC =  α2  x rCO21 +  ω2  x  ω2  x rCO21 =  Matrix([[494.173061379222, -5648.42393818020, 0]])
    AD3 = AC +  α3  x rDC +  ω3  x  ω3  x rDC =  Matrix([[-8.08074589703349*aa3 + 57.0279088954586, 99.6729729954293*aa3 - 5683.86442712301, 0]])
    AD4 =  α4  x rDO41 +  ω4  x  ω4  x rDO41 =  Matrix([[-77.8143817610846*aa4 - 1274.15138480262, 18.5720755743718*aa4 - 5338.51490542404, 0]])

Solve for aa3 and aa4

```
Eqa1 = sp.Eq(AD3[0],AD4[0])
Eqa2 = sp.Eq(AD3[1],AD4[1])
ansa = sp.solve([Eqa1,Eqa2])
```

Print the answer.

```
print("Solution for aa3 and aa4 of matching accelerations equations of point D\n",ansa)
```

    Solution for aa3 and aa4 of matching accelerations equations of point D
     {aa3: 0.282726957185647, aa4: -17.0777511678876}

Substitute the answer in the vectors aa3 and aa4

```
alpha3 = alpha3.subs(aa3,ansa[aa3])
alpha4 = alpha4.subs(aa4,ansa[aa4])
```

```
print("Substituting aa3 and aa4 in vectors",alp_str[1]," and ",alp_str[2],"\n",alp_str[1],"=",alpha3,"\n",alp_str[2],"=",alpha4)
```

    Substituting aa3 and aa4 in vectors α3  and  α4
     α3 = Matrix([[0, 0, 0.282726957185647]])
     α4 = Matrix([[0, 0, -17.0777511678876]])

Calculate the rest of the accelerations.

```
AP3 = AC + alpha3.cross(rPC) + omega3.cross(omega3.cross(rPC))
AF5 = AP3 + alpha5.cross(rFP) + omega5.cross(omega5.cross(rFP))
AF6 = alpha6.cross(rFO61) + omega6.cross(omega6.cross(rFO61))
```

Print the results of the accelerations.

```
print("Point accelerations\nAP3 = AC + ",alp_str[1]," x rPC + ",w_str[1]," x ",w_str[1]," x rPC = ",AP3,"\nAF5 = AP3 + ",alp_str[3]," x rFP + ",w_str[3]," x ",w_str[3]," x rFP = ",AF5,"\nAF6 = ",alp_str[4]," x rFO61 + ",w_str[4]," x ",w_str[4]," x rFO61 = ",AF6)
```

    Point accelerations
    AP3 = AC +  α3  x rPC +  ω3  x  ω3  x rPC =  Matrix([[263.445184231292, -5816.86102785437, 0]])
    AF5 = AP3 +  α5  x rFP +  ω5  x  ω5  x rFP =  Matrix([[15.1079322561889*aa5 - 306.349291009673, 119.045161106793*aa5 - 5744.54883765817, 0]])
    AF6 =  α6  x rFO61 +  ω6  x  ω6  x rFO61 =  Matrix([[-96.2492132913208*aa6 + 1637.83554061634, -27.130959083708*aa6 - 5810.35052238718, 0]])

Solve for aa5 and aa6.

```
Eqa3 = sp.Eq(AF5[0],AF6[0])
Eqa4 = sp.Eq(AF5[1],AF6[1])
ansa2 = sp.solve([Eqa3,Eqa4])
```

Print the answer.

```
print("Solution for aa5 and aa6 of matching accelerations equations of point F\n",ansa2)
```

    Solution for aa5 and aa6 of matching accelerations equations of point F
     {aa5: -5.34760762788408, aa6: 21.0388849545367}

Substitute the answer in the vectors aa5 and aa6

```
alpha5 = alpha5.subs(aa5,ansa2[aa5])
alpha6 = alpha6.subs(aa6,ansa2[aa6])
```

```
print("Substituting aa5 and aa6 in vectors",alp_str[3]," and ",alp_str[4],"\n",alp_str[3],"=",alpha5,"\n",alp_str[4],"=",alpha6)
```

    Substituting aa5 and aa6 in vectors α5  and  α6
     α5 = Matrix([[0, 0, -5.34760762788408]])
     α6 = Matrix([[0, 0, 21.0388849545367]])