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

---

## Solution

---

# Spring mass system.


The second order equation of an unforced spring/mass system with mass m, damping constant b and spring constant k is given by:<p>
$$ mx''(t) + bx'(t) +kx(t) = 0. $$
    
    
Recall that $x(t)$ represents the displacement of the mass from its resting position ( with x>0 the spring is stretched and x<0 the spring is compressed). Therefore the velocity of the moving mass is $ x'(t)$ and the acceleration is $x’’(t)$. To approximate the solution for this system using Euler’s method we convert the equation into a system by introducing the variable
                                                                                                                                
 $$ y = x’.$$
                                                                                                                                
So our two unknowns are: $x$ for the position of the mass and $y$ for the velocity of the mass. As a system the second order equation becomes ( using m = 1)
  
\begin{align}
x'(t) &=  \; y(t) \\
y'(t) &= -k \cdot x(t) - b \cdot y(t)
\end{align}                                                                                                                        

                                                                                                 
and unlike the Predator-Prey model of exercise 2. this system is linear and can be represented in matrix notation:
                                                                                                                                
                                                                                                                                
\begin{align}
        \begin{pmatrix}
        x'(t) \\
        y'(t) 
        \end{pmatrix}=       
        \begin{pmatrix}
        0 & 1 \\
        -k & -b 
        \end{pmatrix}\cdot
        \begin{pmatrix}
        x(t) \\
        y(t) 
        \end{pmatrix}
    \end{align}
                                                                                                                                
                                                                                                                                
   ---

### Collaboration.  Students are allowed to work together on this project. Colaboration is encouraged. 
 However you final submission must be your own work.
  


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


```python
# import libraries

import numpy as np
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.

For this exercise you will approximate the solution of an unforced mass-spring system with
mass m =1, damping constant b = 1 and spring constant k = 2. The second order equation for this system is given by:

$$x''(t) + 1x'(t) + 2x(t) = 0.$$

Use the initial conditions $x(0) = 0$ and $ x'(0)=y(0) = 1$. Use the improved Euler’s method to solve the system. You will need to convert the second order equation to a system as demonstrated above. Choose a Δt=h of 0.1 and a time interval of 12 seconds. Create two plots. One showing x(t) versus time and the second showing the solution curve "trajectory" for the system in the xy (phase) plane.



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


```python
# define the slope functions

def f(x,y):
    return y

def g(x,y):
    return -1*y -2*x


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

# Initialeze arrays to store the results.

x = np.zeros_like(t)
y = np.zeros_like(t)

# Set initial conditions:

x[0] = 0
y[0] = 1

# implement Euler's method

for i in range (len(t)-1):
    
    F1 = f(x[i],y[i])
    G1 = g(x[i],y[i])
    
    F2 = f(x[i]+F1*h,y[i]+G1*h)
    G2 = g(x[i]+F1*h,y[i]+G1*h)
    
    slope1 = (F1+F2)/2
    slope2 = (G1+G2)/2
    
    x[i+1] = x[i] + slope1 * h
    y[i+1] = y[i] + slope2 * h
    
    

```

### Execute the cell below to graph your results.
You should see a damped sinusoid.


```python
plt.figure(figsize=(8,5))
plt.plot(t,x)
plt.grid()
plt.xlabel('t')
plt.ylabel('x');
plt.title("Damping b = 1");
```

### Execute the cell below to see the trajectory in the phase plane for this problem. 
You should see a spiral.


```python
plt.figure(figsize=(8,5))
plt.plot(x,y,linewidth=2)
plt.xlim(-1.5,1.5)
plt.ylim(-1.5,1.5)
plt.grid()
plt.xlabel('x-displacement')
plt.ylabel('y-velocity');
plt.title("Trajectory for spring mass system b = 1");
```

### Repeat the steps above for b=3 and b=0. Copy and paste cells above and make the necessary modifications.  You should show time plots and phase plane plots for b=0 and b=3.



```python
# b = 0

# define the slope functions

def f(x,y):
    return y

def g(x,y):
    return  -2*x


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

# Initialeze arrays to store the results.

x = np.zeros_like(t)
y = np.zeros_like(t)

# Set initial conditions:

x[0] = 0
y[0] = 1

# implement Euler's method

for i in range (len(t)-1):
    
    F1 = f(x[i],y[i])
    G1 = g(x[i],y[i])
    
    F2 = f(x[i]+F1*h,y[i]+G1*h)
    G2 = g(x[i]+F1*h,y[i]+G1*h)
    
    slope1 = (F1+F2)/2
    slope2 = (G1+G2)/2
    
    x[i+1] = x[i] + slope1 * h
    y[i+1] = y[i] + slope2 * h
    
    

```


```python
plt.figure(figsize=(8,5))
plt.plot(t,x)
plt.grid()
plt.xlabel('t')
plt.ylabel('x');
plt.xticks(np.arange(0,13,.5))
plt.title("Damping b = 0");
```


```python
plt.figure(figsize=(8,5))
plt.plot(x,y,linewidth=2)
plt.xlim(-1.5,1.5)
plt.ylim(-1.5,1.5)
plt.grid()
plt.xlabel('x-displacement')
plt.ylabel('y-velocity');
plt.title("Trajectory for spring mass system b = 0");
```


```python
# b=3

# define the slope functions


def f(x,y):
    return y

def g(x,y):
    return -3*y -2*x


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

# Initialeze arrays to store the results.

x = np.zeros_like(t)
y = np.zeros_like(t)

# Set initial conditions:

x[0] = 0
y[0] = 1

# implement Euler's method

for i in range (len(t)-1):
    
    F1 = f(x[i],y[i])
    G1 = g(x[i],y[i])
    
    F2 = f(x[i]+F1*h,y[i]+G1*h)
    G2 = g(x[i]+F1*h,y[i]+G1*h)
    
    slope1 = (F1+F2)/2
    slope2 = (G1+G2)/2
    
    x[i+1] = x[i] + slope1 * h
    y[i+1] = y[i] + slope2 * h
    
    

```


```python
plt.figure(figsize=(8,5))
plt.plot(t,x)
plt.grid()
plt.xlabel('t')
plt.ylabel('x');
plt.title("Damping b = 3");
```


```python
plt.figure(figsize=(8,5))
plt.plot(x,y,linewidth=2)
plt.xlim(-1.5,1.5)
plt.ylim(-1.5,1.5)
plt.grid()
plt.xlabel('x-displacement')
plt.ylabel('y-velocity');
plt.title("Trajectory for spring mass system b = 3");
```

---
# Exercise 2.

For this exercise you will approximate the solution of an undamped periodically forced spring-mass system with mass m =1 and spring constant k = 1, no damping. The second order equation for this system is given by:

$$ x''(t) + x(t) = Cos(\omega t).$$

Use the initial conditions $x(0) = 0$ and $ x'(0) = y(0) = 0$. Use the improved Euler’s method to solve the system. Note this system is not autonomous so the slope functions could depend on x,y and t.  Generate plots of $x(t)$ versus time for ω = 1.1 and ω = 1. Use a time interval of 50 seconds for $ \omega = 1 $ and 150 seconds for $ \omega = 1.1$. For both cases set Δt=0.1. Duplicate and modify the code above to generate the plots.  Label the plots appropriately. There are no phase plane plots for Exercise 2. 

## Solution for $ \omega = 1. $


```python
# define the slope functions

def f(x,y):
    return y

def g(x,y,t):
    return  -1*x + np.cos(t)


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

# Initialeze arrays to store the results.

x = np.zeros_like(t)
y = np.zeros_like(t)

# Set initial conditions:

x[0] = 0
y[0] = 0

# implement Euler's method

for i in range (len(t)-1):
    
    F1 = f(x[i],y[i])
    G1 = g(x[i],y[i],t[i])
    
    F2 = f(x[i]+F1*h,y[i]+G1*h)
    G2 = g(x[i]+F1*h,y[i]+G1*h,t[i+1])
    
    slope1 = (F1+F2)/2
    slope2 = (G1+G2)/2
    
    x[i+1] = x[i] + slope1 * h
    y[i+1] = y[i] + slope2 * h
    
    

```


```python
plt.figure(figsize=(8,5))
plt.plot(t,x)
plt.grid()
plt.xlabel('t')
plt.ylabel('x');
plt.title("Forced Spring-Mass system no damping $\omega = 1.$");
```

## Solution for $ \omega = 1.1 $


```python
# define the slope functions

def f(x,y):
    return y

def g(x,y,t):
    return  -1*x + np.cos(1.1*t)


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

# Initialeze arrays to store the results.

x = np.zeros_like(t)
y = np.zeros_like(t)

# Set initial conditions:

x[0] = 0
y[0] = 0

# implement Euler's method

for i in range (len(t)-1):
    
    F1 = f(x[i],y[i])
    G1 = g(x[i],y[i],t[i])
    
    F2 = f(x[i]+F1*h,y[i]+G1*h)
    G2 = g(x[i]+F1*h,y[i]+G1*h,t[i+1])
    
    slope1 = (F1+F2)/2
    slope2 = (G1+G2)/2
    
    x[i+1] = x[i] + slope1 * h
    y[i+1] = y[i] + slope2 * h
    
    

```


```python
plt.figure(figsize=(12,5))
plt.plot(t,x)
plt.grid()
plt.xlabel('t')
plt.ylabel('x');
plt.title("Forced Srping-Mass system no damping $\omega = 1.1$");
```

## Answer the questions below.  Create a text cell to enter your answers.
1. Expain and discuss the results of exercise 1. 

    For exercise \#1 we see the spring behaviour under three damping conditions. <p>
    b = 0 ~ we see no damping and the spring oscillates.<p>
    b = 1 ~ we see an underdamped behaviour and the oscillation decays.<p>
    b = 3 ~ we see overdamped and the spring returns to equilibrium with no oscillations.<p>
    
2. Compute the period of oscillation for exerice 1. b=0.  Compare to value obtained from the graph.
        
        The period is 2*pi/omega = sqrt(2)*pi = 4.44.  This is close to the value from the graph.
       
        
3. Discuss the results of exercise 2.  Hint: look at the envelopes.
        
        For omega = 1, we are exciting the system at the resonant frequency and thus see RESONANCE, that is the osicllation are increasing. 
        
        For omega = 1.1, the exciting frequencey is "close" to the resonant frequency and we observe the BEATING phenomenon.


```python

```

de4219a648b815a085c85658ea1389dc5c56cdc5
ipynb
Jupyter Notebook
Project#3_solution.ipynb
rmartin977/math---267-Spring-2022
828fce843795318fb1ec32e4dd073b67861e06cf
MIT
Project#3_solution.ipynb
rmartin977/math---267-Spring-2022
828fce843795318fb1ec32e4dd073b67861e06cf
MIT
Project#3_solution.ipynb
rmartin977/math---267-Spring-2022
828fce843795318fb1ec32e4dd073b67861e06cf
MIT
Qwen/Qwen-72B
1. YES
2. YES