<h1 style='font-size:4rem;color:orange;'>Math 267 Project #4

# Dynamical Systems


In simplest terms a Dynamical System is a point moving in space. We call this space a phase space
and the path for the point a trajectory or orbit. The point could represent an actual particle that is moving in physical space or it could represent the instantaneous state of our Dynamical System. The trajectory thus illustrates how the system evolves with time. The evolution follows a vector field. The trajectories could be fixed points, spirals, closed loops that represent periodic behavior or the trajectory could be a chaotic orbit. Both the predator-prey system and the mass-spring system are examples of Dynamical Systems. For the predator-prey system a point the phase space (the xy plane) represents the current population of the foxes and rabbits. For the mass spring system a point represents the current position and velocity of the mass. Thus for this case the instantaneous energy in the system is a function of the point. More on this later..

---

# The Pendulum

The equation of motion for a driven damped pendulum with mass 1 and arm length 1 is:

$$ \theta''(t) + b\; \theta'(t) + g \sin(\theta) = A \sin(\omega_f  t)$$

where $\theta$ is the angular displacement for the pendulum. The forcing term is a sinusoidal driving force with an amplitude of A and a period of $2 \pi / \omega_f$. We derive this equation in class. Please refer to your text for further details. It is assumed that the pendulum is attached to a wall such that it is free to rotate through 360 degrees, thus $\omega + 2\pi$ and $\omega$ correspond to the same pendulum position. We can reduce this second order equation to a system by introducing the variable $\omega = \theta’$ , where $\omega$ represents the angular velocity for the pendulum motion. For convenience let us assume that units are chosen so that g = 1. The system is written as:

\begin{align}
\theta'(t) &=  \; \omega(t) \\
\omega'(t) &= -b \omega(t) + A \sin(\omega_ft)
\end{align} 




## 1. The unforced and undamped Pendulum.

For the case where there is both no damping or external forcing, the equation for the pendulum
becomes:

$$\theta'(t) + \sin(\theta) = 0$$


This equation while simple in appearance cannot be solve analytically. For small $\theta$ we could use the small angle approximation, namely $\sin(\theta) ≈\theta$, and the equation becomes linear. Otherwise numerical methods are required to solve the system. We can however get an idea of what the solutions should look like. With no damping or means to dissipate the energy, the energy E will be constant along all solutions of the system. The total energy in the system, the potential and kinetic energy, is given by the energy function:

$$E(\theta,\omega) = \frac{1}{2} \omega^2 + 1 - \cos(\theta)$$.


Each trajectory in the phase plane is confined to one or another of the level curves for the energy function. So by plotting the contour map for the energy function we get an idea of what the trajectories should look like. A plot of the contour map for the energy function is shown in the figure below. If we add directions to the level curves we have a phase portrait for the pendulum system. Note that only a unique trajectory can pass through each point in the phase plane. Why?




To determine the direction for these trajectories we simply look at the vector field for the system. Writing the unforced, undamped equation as a system gives

\begin{align}
\theta'(t) &=  \; f(\theta, \omega) =\omega \\
\omega'(t) &= \;g(\theta,\omega)= -\sin(\omega).
\end{align} 


So our vector field is:

$$\mathbf{F}(θ,ω) = (f(θ,ω),g(θ,ω)) = (ω,-sin(θ)).$$

Evaluating the vector field at the point (0,1) gives $F(0,1) = (1,0)$ and the trajectory is thus pointing to the right at this point meaning that $\theta$ is increasing. Repeating this for other points in the phase plane allows us to attach direction arrows to the contour curves. The phase portrait for the undamped pendulum is shown in the figure below.




The solution curves encircling the origin where $-π < θ(t) < π$ correspond to the pendulum oscillating about the downward rest position. The solution curve at the top of the phase portrait corresponds to the pendulum spinning forever in the clockwise direction. The opposite is true for the trajectory at the bottom of the phase portrait.

---

## 2. The unforced and damped Pendulum


For the case where there is damping but no external forcing, the equations for the pendulum becomes:

\begin{align}
\theta'(t) &=  \; f(\theta, \omega) =\omega \\
\omega'(t) &= \;g(\theta,\omega)= -\sin(\theta)-b\omega.
\end{align} 

The equilibrium points for this plane autonomous system are determined by setting:

$$\mathbf{F}(θ, ω) = 0 \;or f (θ, ω) = g(θ, ω) = 0$$

Solving gives the points (0, 0) and (0, π).

These points correspond to the pendulum being positioned either in the downward rest position or straight up (θ=π). To determine the stability of these points we need to linearize the system by computing the Jacobian matrix. See 10.3 in the text for further information on the Jacobian. The Jacobian matrix A for the system is defined as a matrix of partial derivatives.



\begin{align} A=             
        \begin{pmatrix}
        \frac{\partial}{\partial \theta}f(\theta,\omega)
         & \frac{\partial}{\partial \omega}f(\theta,\omega)\\
        \frac{\partial}{\partial \theta}g(\theta,\omega)
         & \frac{\partial}{\partial \omega}g(\theta,\omega)
         \end{pmatrix}
        \end{align}
        
        
For the undriven, damped pendulum system the Jacobian is:
        
\begin{align} A=             
        \begin{pmatrix}
        0 & 1\\
        -\cos(\theta)
         & -b
         \end{pmatrix}
        \end{align}
        
Now to determine the stability of an equilibrium point, the Jacobian is evaluated at that point and the eigenvalues are tested. For example, for the equilibrium point (0,0) (the rest position) the linearized system is:

\begin{align}
        \begin{pmatrix}
        \theta'(t) \\
        \omega'(t) 
        \end{pmatrix}=       
        \begin{pmatrix}
        0 & 1 \\
        -1 & -b 
        \end{pmatrix}\cdot
        \begin{pmatrix}
        \theta(t) \\
        \omega(t) 
        \end{pmatrix}
\end{align}

With a damping constant b = 0, the eigenvalues are pure imaginary and the origin is a “center” thus the trajectories will be closed loops with the origin as their center. For b>0 both eigenvalues will have negative real parts and the point (0,0) will be an attractor. Which is what we expect.

---


## 3. The forced and damped Pendulum


For the case where there is damping and external forcing, the equations for the pendulum becomes:

\begin{align}
\theta'(t) &=\omega \\
\omega'(t) &= \; -\sin(\theta)-b\omega + A \sin(\omega_f t).
\end{align} 

This is where things get interesting and we will see that for certain parameter values the pendulum will exhibit chaotic behavior. Sensitive dependence on initial conditions is a hallmark of chaos. The plot below shows the response of a pendulum to the following initial conditions. For s(t) the initial angular displacement θ(0) = 0 and for u(t) the initial displacement is θ(0) = 0.001. We notice the curves diverge considerably after about 40 seconds. This is sensitive dependence on initial conditions.



---

### Run the cell below to import necessary libraries.


```python
# import libraries

import numpy as np
from numpy import sin,pi
import matplotlib.pyplot as plt
from IPython.display import Image

# uncomment the line below if you are running a macbook
#%config InlineBackend.figure_format ='retina'
```

## Exericse 1.a

For this exercise you will approximate the solution of an unforced, undamped defined by the system:

\begin{align}
\theta'(t) &=  \;  \omega \\
\omega'(t) &= \; -\sin(\theta).
\end{align} 

You are to create two plots of $\theta$ versus t. For the first plot solve the sytem with initial condition:

$$
\begin{align}
\theta(0) &=  \;0.1 \\
\omega(0) &= \; 0
\end{align} 
$$

For the second plot solve the system with the initial conditions:

$$
\begin{align}
\theta(0) &=  \;1 \\
\omega(0) &= \; 0
\end{align} 
$$

Use the improved Euler method( Δt=h = 0.1) and a time interval of 12 seconds.


## Complete the code cell below to solve the sytstem.


```python
# define the slope functions

def f(theta,omega,t=None):
    return omega

def g(theta,omega,t=None):
    return -np.sin(theta)


h = 0.1  # set delta t
t = np.arange(0,12,h)

# Initialize arrays to store the results.

theta = np.zeros_like(t)
omega = np.zeros_like(t)

# Set initial conditions:

theta[0] = 0
omega[0] = 0.1

# implement the Improved Euler's method

# Your code goes here, look at project #2
```

### Execute the cell below to graph your results.
You should see a  sinusoid as the pendulum is set into infinite oscillation (no damping).


```python
plt.figure(figsize=(8,5))
plt.plot(t,theta)
plt.grid()
plt.xlabel('t')
plt.ylabel(r'$\theta$');
plt.plot(2*np.pi,0,'o')
plt.title(r"Exercise #1a. $\theta(0)=0\;,\omega(0)=0.1$");
```

### Execute the cell below to graph the trajectory in the phase plane.
You should see a  small circle as the pendulum is set into infinite oscillation (no damping).


```python
plt.figure(figsize=(8,5))
plt.plot(theta,omega,linewidth=2)
plt.xlim(-2*pi,2*pi)
plt.ylim(-4,4)
plt.grid()
plt.xlabel(r'$\theta$-displacement')
plt.ylabel(r'$\omega$');
plt.xticks([-2*pi,-pi,0,pi,2*pi],['-2$\pi$','$-\pi$',0,'$\pi$','$2\pi$'])
plt.title(r"Trajectory for pendulum-Exercise #1a. $\theta(0)=0\;,\omega(0)=0.1$");
```

 ## Exercise 1.b
 
 Repeat the steps above for $\theta(0)=0 \;\text{and}\; \omega(0) = 1.0$ Copy and paste cells above and make the necessary modifications. Make a plot of $\theta$ versus t. Also plot the trajectory in the phase plane. Make the appropriate change to the title of your plot.



```python

```


```python

```

## Exericse 2.

Use the improved Euler’s method to approximate and graph the solution of the above system with the initial conditions θ(0) = 0 and ω(0) =2.5. Create a θ versus t and a phase plot. Your $\theta$ versus t plot will not be a sinusoid. Again plot your results for the interval [0,12]. 



```python

```


```python

```

## Exercise 3.

Use the improved Euler’s method to approximate and graph the solution of the above system with a damping constant of b = 0.5. Use the initial conditions θ(0) = 1 and ω(0) =0. Create a θ versus t plot
and a phase plot θ versus ω. Give each plot the following title: "Exercise #3 damped pendulum b=0.5".


```python
# Your code here.
```


```python

```


```python

```

## Exercise 4

In the following exercises we will approximate and graph the solutions for the forced,damped pendulum system with different damping constants. For each exercise we use  A = 1.5, ωf = 2/3, and initial condition (θ(0) ,ω(0)) = (0,0).

## Exercise 4.a      


Use the improved Euler method to solve the system with b = 0.8, h= Δt=0.1 and time interval of \[0,90]. Make a $\theta$ versus t plot and a phase plot. Give your graphs appropriate titles.

Note: When creating the phase plot we want to omit the transient from the plot and only view the steady state response.  This can be accomplished by plotting just the last third of our data values as follows:

```python

plt.plot(theta[600:900],omega[600:900],linewidth=1)

```

Also note a line width of 1 is recommended for the phase plot.


```python
# define the slope functions

b=0.8

def f(theta,omega,t):
    return omega

def g(theta,omega,t):
    return -np.sin(theta) -b*omega + 1.5*np.sin(2/3*t)

# Your code starts here....


```

### Execute the cell below to graph your results.



```python
plt.figure(figsize=(8,5))
plt.plot(t,theta)
plt.grid()
plt.xlabel('t')
plt.ylabel(r'$\theta$');
plt.title(r"Exercise #4a. b =0.8");
```

### Execute the cell below to see the trajectory in the phase plane for this problem. 



```python
plt.figure(figsize=(10,7))
plt.plot(theta[600:900],omega[600:900],linewidth=1)
plt.xlim(-2*pi,2*pi)
plt.ylim(-4,4)
plt.grid()
plt.xlabel(r'$\theta$-displacement')
plt.ylabel(r'$\omega$');
plt.xticks([-2*pi,-pi,0,pi,2*pi],['-2$\pi$','$-\pi$',0,'$\pi$','$2\pi$'])
plt.title(r"Exercise #4a. b =0.8");
```

## Exercise 4.b

Change b = 0.735 and repeat the steps above. Use h = 0.1 and a time interval of \[0,90].


```python
# Place your code here to solve the system in this cell.
```


```python
# Place your code here to make a theta versus t plot in this cell.
```

### Execute the cell below to see the trajectory in the phase plane for this problem. 



```python
plt.figure(figsize=(10,7))
plt.plot(theta[600:900],omega[600:900],linewidth=1)
plt.xlim(-2*pi,2*pi)
plt.ylim(-4,4)
plt.grid()
plt.xlabel(r'$\theta$-displacement')
plt.ylabel(r'$\omega$');
plt.xticks([-2*pi,-pi,0,pi,2*pi],['-2$\pi$','$-\pi$',0,'$\pi$','$2\pi$'])
plt.title(r"Exercise #4b. b =0.735");
```

## Exercise 5.

Write the necessary code to  recreate the plot that demonstrates sensive dependence on intial condtions. See \#3 above The plot labeled s(t) has an initial value of 0.0 and u(t) has intial value of 0.001. The sinusoid input has an amplitude of 1 and omega = 2/3. The damping for the spring is 0.05. Plot the output for the range (0,80). 

For this exercise you will be solving the system twice.  Just cut and past code above make the necessary modificatons and in the last line of the cell enter the Python 
command.

```python
s=theta;

```
Do the same for a second cell and again end with the command:

```python
u=theta;

```


```python
# Your code here.
```


```python
# Your code here.
```

### Execute the cell below to see a plot of your results



```python
plt.figure(figsize=(8,5))
plt.plot(t,s,label="$s(t)$")
plt.plot(t,u,label="$u(t)$")
plt.grid()
plt.xlabel('t')
plt.legend()
plt.title(r"Exercise #5");
```

## Answer the questions below.  You can edit the text cell to enter your answers.
1. For exercises 1.a and 1.b compute the period of oscillation using the small angle approximation.  Compare this to the value obtained from the graphs and explain the results.

2. Explain the results of exercise 2, that is, explain the shape of the graphs.

3. Look at the phase plots for exercise \#4. Outputs for both b = 0.8 and b = 0.735 are periodic.  How does the period of the output (b=0.735) compare to the output for b=0.8.  

To see your input along with the output add the line to the cell that creates the theta versus t plot.

```python
plt.plot(t,1.5*sin(2/3*t))

```

4. In exercise 5. we see that a very small change in initial conditions produce a dramatic change in output.  Will the same phenomenon occur with the spring-mass system?  In other words can our spring-mass system exhibit sensitive dependence on initial conditions?

03f98f91729e52c9d95b3c3b8aaeefb78879277f
ipynb
Jupyter Notebook
project#4.ipynb
rmartin977/math---267-Spring-2022
828fce843795318fb1ec32e4dd073b67861e06cf
MIT
project#4.ipynb
rmartin977/math---267-Spring-2022
828fce843795318fb1ec32e4dd073b67861e06cf
MIT
project#4.ipynb
rmartin977/math---267-Spring-2022
828fce843795318fb1ec32e4dd073b67861e06cf
MIT
Qwen/Qwen-72B
1. YES
2. YES