Midterm 2 (Due 5 Apr)

Solution Midterm 2, Spring 2024

import numpy as np
from math import *
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
plt.style.use('seaborn-v0_8-colorblind')

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 4
      2 from math import *
      3 import matplotlib.pyplot as plt
----> 4 import pandas as pd
      5 get_ipython().run_line_magic('matplotlib', 'inline')
      6 plt.style.use('seaborn-v0_8-colorblind')

ModuleNotFoundError: No module named 'pandas'

Part 1, The Duffing Oscillator (60pt)

Consider the following equation of motion for a damped and driven oscillator that is not necessarily linear. This is the Duffing equation, which describes a damped driven simple harmonic oscillator with a cubic non-linearity. It’s equation of motion is given by:

\[\ddot{x} + \delta \dot{x} + \alpha x + \beta x^3 = \gamma \cos(\omega t).\]

where \(\alpha\), \(\beta\), \(\delta\), \(\gamma\), and \(\omega\) are constants.

Below is a figure showing the strange attractor of the Duffing Oscillator over four periods. It is locally chaotic, but globally it is stable. This is a common feature of non-linear systems.

We focus on this one dimensional case, but we will analyze it in parts to put together a full picture of the dynamics of the system. If you are looking for useful parameter choices and initial conditions, we suggest using those listed on the Wikipedia page.

Recall that the potential for this system can give rise to different kinds of behavior depending on the parameters. The figure below shows the potential for different choices of \(\alpha\) and \(\beta\). Make sure you choose your parameters to match the potential you want to study - pick a stable well for solving for trajectories.

Undamped and undriven case (20 points)

1. (5 pt) Consider the undriven and undamped case, \(\delta = \gamma = 0\). What potential would give rise to this kind of equation of motion? Demonstrate that potential gives the expected equation of motion.

Solution

We are focusing on the one-dimensional and undamped case, so we can write the equation of motion as:

\[\ddot{x} + \alpha x + \beta x^3 = 0.\]

This equation of motion is related to the spatial dependence of a potential energy. We can see this because the two terms that are not the acceleration are dependent only on position. We can write this as:

\[\ddot{x} = \dfrac{F_{x}}{m} = -\dfrac{1}{m}\dfrac{dV}{dx} = -\alpha x - \beta x^3.\]

This means that the potential energy is, up to a constant, \(V(x) = \dfrac{1}{2} m\alpha x^2 + \dfrac{1}{4} m\beta x^4\).

We can verify that this potential gives the expected equation of motion by taking the derivative of the potential with respect to \(x\) and plugging it into the equation of motion. We find:

\[\dfrac{dV}{dx} = m\alpha x + m\beta x^3,\]

So that the equation of motion becomes:

\[\ddot{x} = -\dfrac{1}{m}\dfrac{dV}{dx} = -\alpha x - \beta x^3.\]

2. (5 pt) Find the equilibrium points of the system for the undriven and undamped case. What are the conditions for stable and unstable equilibrium points? Plot the potential for some choices of \(\alpha\) and \(\beta\) to illustrate your findings.

Solution

Here we can find the equilbrium points by setting the net force to zero (i.e., the derivative of the potential energy is zero). We have:

\[\dfrac{dV}{dx} = m\alpha x + m\beta x^3 = 0.\]

We can factor this to find the equilibrium points:

\[x(\alpha + \beta x^2) = 0.\]

So that there are in principle up to three equilibrium points.

\[x = 0, \quad x = \pm \sqrt{-\dfrac{\alpha}{\beta}}.\]

The stability of these points is determined by the concavity of the potential. If the potential is concave up, the equilibrium point is stable. If the potential is concave down, the equilibrium point is unstable. We can see this by looking at the second derivative of the potential:

\[\dfrac{d^2V}{dx^2} = m\alpha + 3m\beta x^2.\]

For \(x=0\), we have

\[\dfrac{d^2V}{dx^2} = m\alpha,\]

and for \(x = \pm \sqrt{-\dfrac{\alpha}{\beta}}\), we have

\[\dfrac{d^2V}{dx^2} = m\alpha - 3m\beta\dfrac{\alpha}{\beta} = m\alpha - 3m\alpha = -2m\alpha.\]

Note that there’s a condition on the existence on two of the equilibrium points, namely that \(\alpha < 0\) and \(\beta > 0\) or \(alpha > 0\) and \(\beta < 0\). This is because the square root of a negative number is imaginary, and we need a real value for \(x\).

• If \(\alpha > 0\), the equilibrium point at \(x = 0\) is stable.

• If \(\alpha > 0\) and \(\beta < 0\), the equilibrium points at \(x = \pm \sqrt{-\dfrac{\alpha}{\beta}}\) exist and are unstable.

• If \(\alpha < 0\), the equilibrium point at \(x = 0\) is unstable.

• If \(\alpha < 0\) and \(\beta > 0\), the equilibrium points at \(x = \pm \sqrt{-\dfrac{\alpha}{\beta}}\) exist and are stable.

The picture above shows the potential for different choices of \(\alpha\) and \(\beta\). We reproduce that figure below.

## Define the static Duffing Potential
def V(x, alpha, beta):
    return 0.5*alpha*x**2 + 0.25*beta*x**4

## Define a series of parameters for the potential

alpha = [1.0, 1.0, -1.0, -1.0]
beta = [1.0, -1.0, 1.0, -1.0]

## Create four plots of the potential

x = np.linspace(-2,2,1000)
plt.figure(figsize=(12,8))
for i in range(4):
    plt.subplot(2,2,i+1)
    plt.plot(x,V(x,alpha[i],beta[i]))
    plt.xlabel('x')
    plt.ylabel('V(x)')
    plt.title(r'$\alpha = $' + str(alpha[i]) + r', $\beta = $' + str(beta[i]))
    plt.grid()
    
plt.tight_layout()

3. (5 pt) Plot the phase space for the undriven and undamped case. What are the trajectories of the system? What is the behavior of the system? Does it match your expectations from the potential? The figure above only starts this discussion.

Solution

For this we will choose the double well stable potential where \(\alpha < 0\) and \(\beta > 0\). We can plot the phase space by plotting the velocity as a function of position. And we do that for the choices of \(\alpha = -1\) and \(\beta = 1\). We plot that below (with \(m =1\)).

We can see that the system always exhibits oscillatory behavior. The system will somtimes oscillate between the two stable equilibrium points. For lower energy choices, it might oscillate around one or the other stable points. This is consistent with the potential we chose, which has two stable equilibrium points.

def DuffingStatic(X,V,alpha=-1.0,beta=1.0):
    dX = V
    dV = -alpha*X - beta*X**3
    return dX, dV

def DuffingPhaseSpace(x_lim, y_lim, grid_size):
    
    x = np.linspace(-x_lim, x_lim, grid_size)
    y = np.linspace(-y_lim, y_lim, grid_size)
    
    X, Y = np.meshgrid(x, y)
    
    dX, dY = DuffingStatic(X, Y)
    
    return X, Y, dX, dY

# Generate the phase space for the Duffing Oscillator
x_lim = 2
y_lim = 2
grid_size = 20

X, Y, dX, dY = DuffingPhaseSpace(x_lim, y_lim, grid_size)

fig = plt.figure(figsize=(12,8))

plt.streamplot(X, Y, dX, dY, color='b', linewidth=1, density=2)
plt.xlabel('x')
plt.ylabel('v')
plt.title('Duffing Oscillator Phase Space')
plt.grid()

4. (5 pt) Here you are likely to need a numerical solver. Solve the equation of motion for the undriven and undamped case for a reasonable choice of \(\alpha\) and \(\beta\) and initial conditions. Plot the position as a function of time. What is the behavior of the system? Does it match your expectations from the potential and the phase portrait?

Solution

We will solve the equation of motion for two cases, one where the oscillator moves between the two stable points and one where it stays near one or the other. Here it’s a measure of the total energy (which is conserved) that determines the behavior of the system. We will continue to use \(\alpha = -1\) and \(\beta = 1\).

But we generate three scenarios. One where the total energy is high, one where the total energy is low, and one where the total energy is in between. We plot the position as a function of time for these three cases below.

• For high energy, the system oscillates between the two stable points. We choose initial conditions \(x(0) = 1\) and \(\dot{x}(0) = 1.5\).

• For low energy, the system oscillates around one of the stable points. We choose initial conditions \(x(0) = 1.0\) and \(\dot{x}(0) = 0.3\). And we choose additional initial conditions \(x(0) = -1.0\) and \(\dot{x}(0) = 0.3\), to see the system oscillate around the other stable point.

• For intermediate energy, the system oscillates around one of the stable points. We choose initial conditions \(x(0) = 0.0\) and \(\dot{x}(0) = 0.01\).

We write the code below.

These different trajectories are consistent with the potential we chose, which has two stable equilibrium points, and we can see that the system oscillates between the two stable points for high energy, and around one (or the other) of the stable points for low energy.

## Euler Cromer First Order Solver for the Duffing Oscillator

def DuffingEulerCromer(x0, v0, alpha=-1.0, beta=1.0, dt=0.001, t_max=30):
    
    t = np.arange(0, t_max, dt)
    
    x = np.zeros_like(t)
    v = np.zeros_like(t)
    
    x[0] = x0
    v[0] = v0
    
    for i in range(1, len(t)):
        x[i] = x[i-1] + v[i-1]*dt
        v[i] = v[i-1] - (alpha*x[i-1] + beta*x[i-1]**3)*dt
        
    return t, x, v

## Solve the Duffing Oscillator for our set of initial conditions

## Define the initial conditions
initialPositions = [1, 1.0, -1.0, 0.0]
initialVelocities = [1.5, 0.3, 0.3, 0.01]

## Find numerical solutions for the Duffing Oscillator
## for our choices of initial conditions
## store the results in separate pandas dataframes

t, x, v = DuffingEulerCromer(initialPositions[0], initialVelocities[0])
highEnergyDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

t, x, v = DuffingEulerCromer(initialPositions[1], initialVelocities[1])
lowEnergy1DF = pd.DataFrame({'t': t, 'x': x, 'v': v})

t, x, v = DuffingEulerCromer(initialPositions[2], initialVelocities[2])
lowEnergy2DF = pd.DataFrame({'t': t, 'x': x, 'v': v})

t, x, v = DuffingEulerCromer(initialPositions[3], initialVelocities[3])
medEnergyDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

## Plot the numerical solutions for the Duffing Oscillator in time
fig = plt.figure(figsize=(12,8))

plt.plot(highEnergyDF['t'], highEnergyDF['x'], label='High Energy')
plt.plot(lowEnergy1DF['t'], lowEnergy1DF['x'], label='Low Energy 1')
plt.plot(lowEnergy2DF['t'], lowEnergy2DF['x'], label='Low Energy 2')
plt.plot(medEnergyDF['t'], medEnergyDF['x'], label='Medium Energy')

plt.xlabel('t')
plt.ylabel('x')
plt.title('Duffing Oscillator Solutions')
plt.grid()
plt.legend()

## Plot the numerical solutions for the Duffing Oscillator in phase space

fig = plt.figure(figsize=(12,8))

plt.plot(highEnergyDF['x'], highEnergyDF['v'], label='High Energy')
plt.plot(lowEnergy1DF['x'], lowEnergy1DF['v'], label='Low Energy 1')
plt.plot(lowEnergy2DF['x'], lowEnergy2DF['v'], label='Low Energy 2')
plt.plot(medEnergyDF['x'], medEnergyDF['v'], label='Medium Energy')

plt.xlabel('x')
plt.ylabel('v')
plt.title('Duffing Oscillator Phase Space Trajectories')
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x14103dd90>

Damped and undriven case (20 points)

5. (5 pt) Now consider the damped and undriven case, \(\gamma = 0\). What is the equation of motion in this case? Can you develop this equation of motion fully from a potential (like we did above)? If so, what is the potential? If not, why not?

Solution

In this case, the equation of motion is:

\[\ddot{x} + \delta \dot{x} + \alpha x + \beta x^3 = 0.\]

We are unable to write this full equation of motion as a derivative of a potential energy. This is because the damping term is a velocity dependent term, and the potential energy is only dependent on position. This system is not conservative. We are stuck with the equation of motion as it is.

6. (5 pt) Plot the phase space for the damped and undriven case. What kinds of trajectories are there?

Solution

We can follow the same idea as we did for the static case and plot the phase space. This is possible, because there is no explicit time dependence in the equation of motion. We can plot the phase space for \(\alpha = -1\), \(\beta = 1\), and \(\delta = 0.5\). We plot that below.

We see from that phase space that we should expect every trajectory to spiral inwards towards one or the other stable equilibrium point. This seems to make sense as we are shedding energy due to the damping term, so the system eventually comes to rest at one of the stable points.

def DuffingDamped(X, V, alpha=-1.0, beta=1.0, gamma=0.5):
    dX = V
    dV = -alpha*X - beta*X**3 - gamma*V
    return dX, dV


def DuffingDampedPhaseSpace(x_lim, y_lim, grid_size):
    
    x = np.linspace(-x_lim, x_lim, grid_size)
    y = np.linspace(-y_lim, y_lim, grid_size)
    
    X, Y = np.meshgrid(x, y)
    
    dX, dY = DuffingDamped(X, Y)
    
    return X, Y, dX, dY

# Generate the phase space for the Duffing Oscillator
x_lim = 2
y_lim = 2
grid_size = 20

X, Y, dX, dY = DuffingDampedPhaseSpace(x_lim, y_lim, grid_size)

fig = plt.figure(figsize=(12,8))

plt.streamplot(X, Y, dX, dY, color='b', linewidth=1, density=2)
plt.xlabel('x')
plt.ylabel('v')
plt.title('Duffing Oscillator Phase Space')
plt.grid()

7. (10 pt) Here you are likely to need a numerical solver. Solve the equation of motion for the damped and undriven case for a reasonable choice of \(\alpha\), \(\beta\), \(\delta\), and initial conditions. Plot the position as a function of time. What is the behavior of the system? Does it match your expectations from the phase portrait?

Solution

We will solve the equation of motion for two cases, one where the system spirals inwards towards one of the stable points and one where the system spirals inwards towards the other stable point. We will use \(\alpha = -1\), \(\beta = 1\), and \(\delta = 0.5\).

• For the left hand stable point, we choose initial conditions \(x(0) = 0.5\) and \(\dot{x}(0) = -1.5\).

• For the right hand stable point, we choose initial conditions \(x(0) = -0.5\) and \(\dot{x}(0) = +1.5\). We use Euler Cromer to integrate the equations of motion and plot the position as a function of time (and their phase trajectories) for these two cases below.

## Euler Cromer First Order Solver for the Duffing Oscillator with damping

def DuffingDampedEulerCromer(x0, v0, alpha=-1.0, beta=1.0, gamma=0.5, dt=0.001, t_max=30):
    
    t = np.arange(0, t_max, dt)
    
    x = np.zeros_like(t)
    v = np.zeros_like(t)
    
    x[0] = x0
    v[0] = v0
    
    for i in range(1, len(t)):
        x[i] = x[i-1] + v[i-1]*dt
        v[i] = v[i-1] - (alpha*x[i-1] + beta*x[i-1]**3 + gamma*v[i-1])*dt
        
    return t, x, v

## Solve the Duffing Oscillator for our set of initial conditions

## Define the initial conditions
initialPositions = [0.5, -0.5]
initialVelocities = [-1.5, 1.5]

## Find numerical solutions for the Duffing Oscillator
## for our choices of initial conditions
## store the results in separate pandas dataframes

t, x, v = DuffingDampedEulerCromer(initialPositions[0], initialVelocities[0])
leftDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

t, x, v = DuffingDampedEulerCromer(initialPositions[1], initialVelocities[1])
rightDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

## Plot the numerical solutions for the Duffing Oscillator in time

fig = plt.figure(figsize=(12,8))

plt.plot(leftDF['t'], leftDF['x'], label='Left')
plt.plot(rightDF['t'], rightDF['x'], label='Right')

plt.xlabel('t')
plt.ylabel('x')
plt.title('Duffing Oscillator Solutions')
plt.grid()
plt.legend()

## Plot the numerical solutions for the Duffing Oscillator in phase space

fig = plt.figure(figsize=(12,8))

plt.plot(leftDF['x'], leftDF['v'], label='Left')
plt.plot(rightDF['x'], rightDF['v'], label='Right')

plt.xlabel('x')
plt.ylabel('v')
plt.title('Duffing Oscillator Phase Space Trajectories')
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x141475b20>

Driven and damped case (20 points)

8. (5 pt) Now consider the damped and driven case, \(\gamma \neq 0\). What is the equation of motion in this case? Can you develop this equation of motion fully from a potential (like we did above)? If so, what is the potential? If not, why not?

Solution

Again, we are unable to write this full equation of motion as a derivative of a potential energy. This is because the driving term is time dependent, and the potential energy is only dependent on position. This system is not conservative. We are stuck with the equation of motion as it is.

We can write it as:

\[\ddot{x} + \delta \dot{x} + \alpha x + \beta x^3 = \gamma \cos(\omega t).\]

9. (5 pt) Is it possible to plot the phase space for the driven and damped case? Why or why not? What solutions are available to you to understand the behavior of the system in phase space? Here we are not looking for you to solve the problem, but to conduct research into how people make sense of the behavior of driven and damped systems.

Solution

It is not possible to plot the phase space for the driven and damped case. This is because the driving term is time dependent, and the phase space is a plot of position versus velocity, but it assumes the phase space is time independent. This is where the Poincaré section comes in. The Poincaré section is a plot of position versus velocity at a fixed time interval. This allows us to see the behavior of the system in phase space, even though the system is time dependent.

10. (10 pt) Here you are likely to need a numerical solver. Solve the equation of motion for the driven and damped case for a reasonable choice of \(\alpha\), \(\beta\), \(\delta\), \(\gamma\), \(\omega\), and initial conditions. Plot the position as a function of time. What is the behavior of the system?

Solution

We will solve the equation of motion for the driven and damped case. We will use \(\alpha = -1\), \(\beta = 1\), \(\delta = 0.5\), \(\gamma = 1.0\), and \(\omega = 2\pi\). We will choose the same initial conditions as before for the initially observed left and right stable points. We use Euler Cromer to integrate the equations of motion and plot the position as a function of time for this case below.

The behavior is very complicated, but we can see that the motion is bounded. This is because the system forced and damped. The combination we chose keeps the system from going off to infinity. The system with our chosen parameters appears chaotic, but it seems globally stable (i.e., bounded to a region of space).

## Euler Cromer First Order Solver for the Duffing Oscillator with damping and oscillating force
def DuffingDampedForcedEulerCromer(x0, v0, alpha=-1.0, beta=1.0, delta=0.5, gamma=1.0, omega=1.0, dt=0.001, t_max=10):
    t = np.arange(0, t_max, dt)
    
    x = np.zeros_like(t)
    v = np.zeros_like(t)
    
    x[0] = x0
    v[0] = v0
    
    for i in range(1, len(t)):
        x[i] = x[i-1] + v[i-1]*dt
        v[i] = v[i-1] - (alpha*x[i-1] + beta*x[i-1]**3 + delta*v[i-1] + gamma*cos(omega*t[i-1]))*dt
        
    return t, x, v

## Solve the Duffing Oscillator for our set of initial conditions

## Define the initial conditions
initialPositions = [0.5, -0.5]
initialVelocities = [-1.5, 1.5]
omega = 2*np.pi
T = 2*np.pi/omega
tmax = 30*T
dt = T/10000

## Find numerical solutions for the Duffing Oscillator
## for our choices of initial conditions
## store the results in separate pandas dataframes

t, x, v = DuffingDampedForcedEulerCromer(initialPositions[0], initialVelocities[0], omega=omega, dt=dt, t_max=tmax)
leftDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

t, x, v = DuffingDampedForcedEulerCromer(initialPositions[1], initialVelocities[1], omega=omega, dt=dt, t_max=tmax)
rightDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

## Plot the numerical solutions for the Duffing Oscillator in time

fig = plt.figure(figsize=(12,8))

plt.plot(leftDF['t'], leftDF['x'], label='Left')
plt.plot(rightDF['t'], rightDF['x'], label='Right')

plt.xlabel('t')
plt.ylabel('x')
plt.title('Duffing Oscillator Solutions')
plt.grid()
plt.legend()

## Plot the numerical solutions for the Duffing Oscillator in phase space

fig = plt.figure(figsize=(12,8))

plt.plot(leftDF['x'], leftDF['v'], label='Left')
plt.plot(rightDF['x'], rightDF['v'], label='Right')

plt.xlabel('x')
plt.ylabel('v')
plt.title('Duffing Oscillator Phase Space Trajectories')
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x1415a2d30>

Extra credit (up to 30 points)

11. (10 pts) The Duffing oscillator can be either a soft or hard spring by changing the sign of the \(\alpha\) term. What happens to the behavior of the system when you change the sign of \(\alpha\)? Can you explain this behavior in terms of the potential?

Solution

In all of our choices, we selected \(\alpha < 0\). This means that the potential is a double well potential. If we change the sign of \(\alpha\), the potential becomes a single well potential. This means that the system will always oscillate around the single stable equilibrium point. We can see this in the phase space plot we made below. It looks similar to the regular damped simple harmonic oscillator.

Let’s drive it as before with \(\gamma = 0.5\) and \(\omega = 1\). We will use \(\alpha = 1\), \(\beta = 1\), \(\delta = 0.5\), and the same initial conditions as before. We use Euler Cromer to integrate the equations of motion and plot the position as a function of time for this case below. And it really does look like a damped simple harmonic oscillator.

def DuffingDampedPhaseSpace(x_lim, y_lim, grid_size):
    
    x = np.linspace(-x_lim, x_lim, grid_size)
    y = np.linspace(-y_lim, y_lim, grid_size)
    
    X, Y = np.meshgrid(x, y)
    
    dX, dY = DuffingDamped(X, Y, alpha=1.0)
    
    return X, Y, dX, dY

# Generate the phase space for the Duffing Oscillator
x_lim = 2
y_lim = 2
grid_size = 20

X, Y, dX, dY = DuffingDampedPhaseSpace(x_lim, y_lim, grid_size)

fig = plt.figure(figsize=(12,8))

plt.streamplot(X, Y, dX, dY, color='b', linewidth=1, density=2)
plt.xlabel('x')
plt.ylabel('v')
plt.title('Duffing Oscillator Phase Space')
plt.grid()

## Solve the Duffing Oscillator for our set of initial conditions

## Define the initial conditions
initialPositions = [0.5, -0.5]
initialVelocities = [-1.5, 1.5]
omega = 2*np.pi
T = 2*np.pi/omega
tmax = 30*T
dt = T/10000

## Find numerical solutions for the Duffing Oscillator
## for our choices of initial conditions
## store the results in separate pandas dataframes

t, x, v = DuffingDampedForcedEulerCromer(initialPositions[0], initialVelocities[0], alpha=1.0, omega=omega, dt=dt, t_max=tmax)
leftDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

t, x, v = DuffingDampedForcedEulerCromer(initialPositions[1], initialVelocities[1], alpha = 1.0, omega=omega, dt=dt, t_max=tmax)
rightDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

## Plot the numerical solutions for the Duffing Oscillator in time

fig = plt.figure(figsize=(12,8))

plt.plot(leftDF['t'], leftDF['x'], label='Left')
plt.plot(rightDF['t'], rightDF['x'], label='Right')

plt.xlabel('t')
plt.ylabel('x')
plt.title('Duffing Oscillator Solutions')
plt.grid()
plt.legend()

## Plot the numerical solutions for the Duffing Oscillator in phase space

fig = plt.figure(figsize=(12,8))

plt.plot(leftDF['x'], leftDF['v'], label='Left')
plt.plot(rightDF['x'], rightDF['v'], label='Right')

plt.xlabel('x')
plt.ylabel('v')
plt.title('Duffing Oscillator Phase Space Trajectories')
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x14192d3a0>

12. (20 pts) A critical tool for understanding the behavior of the Duffing oscillator is the Poincaré section. Can you implement a Poincaré section for the driven and damped case? What does it tell you about the behavior of the system?

Solution

In an analysis where we want to use Poincare sections, we need to keep track of the location and velocity at a fixed time interval. We can do this by integrating the equations of motion and then plotting the position and velocity at a fixed time interval. This allows us to see the behavior of the system in phase space, even though the system is time dependent. It’s common to choose the period of the driving force as the time interval. We set that up below.

from matplotlib.cm import viridis

## Solve the Duffing Oscillator for our set of initial conditions

## Define the initial conditions
initialPos = 0.5
initialVel = -1.5
omega = 2*np.pi
T = 2*np.pi/omega
tmax = 100*T
dt = T/10000

## Find numerical solutions for the Duffing Oscillator
## for our choices of initial conditions
## store the results in separate pandas dataframes

t, x, v = DuffingDampedForcedEulerCromer(initialPos, initialVel, alpha=1.0, omega=omega, dt=dt, t_max=tmax)
trajDF = pd.DataFrame({'t': t, 'x': x, 'v': v})

## Extract the position and velocity data for the left and right solutions
## at regular time intervals associated with the period of the driving force
## With omega = 2*pi, the period is 1, we want to sample every 1/T.

trajDF_sampled = trajDF.iloc[::int(T/dt)]

## Plot the numerical solutions for the Duffing Oscillator in phase space
## with the solutions sampled at regular intervals

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

# Main plot
colors = viridis(np.linspace(0, 1, len(trajDF_sampled)))
ax.scatter(trajDF_sampled['x'], trajDF_sampled['v'], c=colors, marker='o', cmap='viridis')
ax.set_xlabel('x')
ax.set_ylabel('v')
ax.grid(True)

plt.tight_layout()

/var/folders/rh/mck_zfls1nz00lkh8jq3cpv80000gn/T/ipykernel_5907/3734017154.py:33: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  ax.scatter(trajDF_sampled['x'], trajDF_sampled['v'], c=colors, marker='o', cmap='viridis')

Part 2, Signal Analysis (40pt)

We learned that any periodic function can be written using a Fourier series. The expansion we wrote is:

\[f(t) = \dfrac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos(n\omega t) + b_n \sin(n\omega t) \right).\]

where \(\omega = 2\pi/T\) and \(T\) is the period of the function. The coefficients \(a_n\) and \(b_n\) are given by:

\[a_n = \dfrac{2}{T} \int_{-T/2}^{T/2} f(t) \cos(n\omega t) dt, \quad b_n = \dfrac{2}{T} \int_{-T/2}^{T/2} f(t) \sin(n\omega t) dt.\]

For this part of the midterm you will analyze a triangular wave both analytically and numerically.

Consider the following triangular wave:

\[\begin{split}f(t) = \begin{cases} 1 - 2t & 0 \leq t < 1/2 \\ 2t - 1 & 1/2 \leq t < 1 \end{cases}\end{split}\]

which repeats every \(T = 1\).

1. (5 pt) Plot the triangular wave. What is the frequency of the wave? Show this on the plot.

Solution

We can prepare a simple plot of the triangular wave. We plot it below. We can see that the wave repeats every \(T = 1\).

T = 1
t_array = np.arange(0, 5*T, T/100)

def f(t, T):
    if t % T < 0.5*T:
        return 1-2*(t % T)/T
    else:
        return 2*(t % T)/T-1
    
wave = np.zeros(len(t_array))

i = 0

for t in t_array:
    
    wave[i]= f(t, T)
    #print('wave:', i, wave[i])
    i += 1
    
plt.figure(figsize=(12,8))

plt.plot(t_array, wave)

plt.xlabel('t')
plt.ylabel('f(t)')

plt.title('Triangle Wave Function')
plt.grid()

2. (15 pt) Determine the Fourier coefficients \(a_n\) and \(b_n\) for the triangular wave. This should be done analytically.

Solution

The original function is defined piecewise,

\[\begin{split}f(t) = \begin{cases} 1 - 2t & 0 \leq t < 1/2 \\ 2t - 1 & 1/2 \leq t < 1 \end{cases}\end{split}\]

so we will eventually have to perform two integrals to find the Fourier coefficients. We can write the Fourier series as:

\[f(t) = \dfrac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos(n\omega t) + b_n \sin(n\omega t) \right).\]

where \(\omega = 2\pi/T = 2\pi\). The coefficients \(a_n\) and \(b_n\) are given by:

\[a_n = \dfrac{2}{T} \int_{T} f(t) \cos(n\omega t) dt, \quad b_n = \dfrac{2}{T} \int_{T} f(t) \sin(n\omega t) dt.\]

We should note that this signal appears to be an even functions, so that we expect only the cosine terms to be present. So we start with the \(a_n\) terms,

\[a_n = \dfrac{2}{1} \int_{0}^{1} f(t) \cos(n2\pi t) dt,\]

\[a_n = I_1 + I_2 = 2 \int_{0}^{1/2} (1 - 2t) \cos(2n\pi t) dt + 2 \int_{1/2}^{1} (2t - 1) \cos(2n\pi t) dt.\]

We evaluate each integral separately. For the first integral, \(I_1\), we have

\[I_1 = 2 \int_{0}^{1/2} (1 - 2t) \cos(2n\pi t) dt = 2 \left[\int_{0}^{1/2} \cos(2n\pi t) dt\right] - 4 \left[\int_{0}^{1/2} t \cos(2n\pi t) dt\right],\]

\[I_1 = \left(\dfrac{\sin(n\pi)}{n\pi}\right) - \left(\dfrac{n\pi \sin(n \pi) + \cos(n\pi) -1}{n^2\pi^2}\right)\]

We should note that \(\sin(n\pi) = 0\) and \(\cos(n\pi) = (-1)^n\) for all integers \(n\). So we have,

\[I_1 = \dfrac{1 - (-1)^n}{n^2\pi^2}.\]

For the second integral, \(I_2\), we have

\[I_2 = 2 \int_{1/2}^{1} (2t - 1) \cos(2n\pi t) dt = 2 \left[\int_{1/2}^{1} 2t \cos(2n\pi t) dt\right] - 2 \left[\int_{1/2}^{1} \cos(2n\pi t) dt\right],\]

\[I_2 = 2 \left(\dfrac{2\sin(n\pi) - 2n\pi \cos(n\pi) - 1}{4n^2\pi^2}\right) - 2 \left(\dfrac{\sin(2n\pi) - \sin(n\pi)}{2n\pi}\right),\]

\[I_2 = \dfrac{1 - (-1)^n}{n^2\pi^2}.\]

So that the \(a_n\) terms are

\[a_n = \left(\dfrac{1 - (-1)^n}{n^2\pi^2} + \dfrac{1 - (-1)^n}{n^2\pi^2}\right) = 2\dfrac{1 - (-1)^n}{n^2\pi^2}.\]

As we indicated above, the \(b_n\) terms are likely zero. We can verify this by calculating them. We have

\[b_n = \dfrac{2}{1} \int_{0}^{1} f(t) \sin(n2\pi t) dt,\]

\[b_n = 2 \int_{0}^{1/2} (1 - 2t) \sin(2n\pi t) dt + 2 \int_{1/2}^{1} (2t - 1) \sin(2n\pi t) dt.\]

We evaluate each integral separately and quickly. For the first integral, we have

\[2 \int_{0}^{1/2} (1 - 2t) \sin(2n\pi t) dt = -\dfrac{\sin(n\pi)-n\pi}{n^2\pi^2} = \dfrac{1}{n\pi}.\]

For the second integral, we have

\[2 \int_{1/2}^{1} (2t - 1) \sin(2n\pi t) dt = -\dfrac{\sin(n \pi) - \sin(2n\pi)+ n\pi\cos(2n\pi)}{n^2\pi^2} = -\dfrac{1}{n\pi}.\]

So that the \(b_n\) terms are zero.

\[b_n = \dfrac{1}{n\pi} - \dfrac{1}{n\pi} = 0.\]

3. (10 pt) What is the expression for the approximation of the triangular wave using the first few terms of the Fourier series?

Solution

The Fourier series for this triangular wave is

\[f(t) = \dfrac{a_0}{2} + \sum_{n=1}^{\infty} a_n \cos(n2\pi t),\]

where \(a_n = 2\dfrac{1 - (-1)^n}{n^2\pi^2}\). The approximation of the triangular wave using the first few terms of the Fourier series is

\[f(t) = \dfrac{1}{2} - \dfrac{2}{\pi^2} \sum_{n=1}^{\infty} \dfrac{1 - (-1)^n}{n^2} \cos(2n\pi t).\]

4. (10 pt) Plot the Fourier series for the first few terms. How many terms do you need to get a good approximation of the triangular wave? Do you notice any problems with the Fourier series? If so, what are they?

Solution

We wrote a few small functions to compute the coefficients below and plot the approximation. From that we notice that there’s maybe a small issue where the triangle wave is not perfectly continuous at the transition points. This is a common issue with Fourier series. We can see that the approximation is good, but not perfect. We zoom in on that part in the second plot.

## Function to calculate the Fourier Series Coefficients
def a0_coeff():
    return 0.5;

def an_coeff(n):
    return 2/(n*np.pi)**2 * (1-(-1)**n)

def bn_coeff(n):
    return 0

## Function to calculate the Fourier Series
def fourier_series(t, T, n_max):
    f = 0
    for n in range(1, n_max+1):
        f += an_coeff(n)*np.cos(2*np.pi*n*t/T) + bn_coeff(n)*np.sin(2*np.pi*n*t/T)
    return f+a0_coeff()

T = 1
n_terms = [1, 2, 3, 4, 5, 10, 20, 50, 100]
t_array = np.arange(0, 5*T, T/100)

plt.figure(figsize=(12,8))
plt.plot(t_array, wave, label='Original Wave')

for n in n_terms:
    f = fourier_series(t_array, T, n)
    plt.plot(t_array, f, label='n = ' + str(n))

plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Fourier Series Approximation of Triangle Wave')
plt.xlim(0,2)
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x1412ddc40>

T = 1
n_terms = [1, 2, 3, 4, 5, 10, 20, 50, 100]
t_array = np.arange(0, 5*T, T/100)

plt.figure(figsize=(12,8))
plt.plot(t_array, wave, label='Original Wave')

for n in n_terms:
    f = fourier_series(t_array, T, n)
    plt.plot(t_array, f, label='n = ' + str(n))

plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Fourier Series Approximation of Triangle Wave (zoomed in)')
plt.xlim(0.95,1.05)
plt.ylim(0.95,1.01)
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x1419e1a00>