# Lecture 3

In this lecture we move onto non-homogeneous second-order ordinary differentials. We want to solve equations of the form

$$
a \frac{d^{2}y}{dx^{2}} + b \frac{dy}{dx} + c y = f(x)
$$

where $a$, $b$ and $c$ are real constants. We've seem previously that for the case $b^{2} \ne 4ac$ and
$f(x) = 0$, the general solution to the above equation is 

$$
y = A_{1} e^{\lambda_{1}} + A_{2}e^{\lambda_{2}}
$$

where $\lambda$ are the roots to the characteristic equation $a\lambda^{2} + b \lambda + c = 0$. For the case that $b^{2} = 4ac$ (repeated roots), the general solution is of the form

$$
y =(A_{1} + A_{2} x) e^{\lambda x}
$$

Recall that the constants $A_{1}$ and $A_{2}$ are determined from the boundary conditions.

### Initialising SymPy

To solve equations using SymPy, we again need to import SymPy, and we;ll call `init_printing()` to get nicely typeset equations:


```
from sympy import *

# This initialises pretty printing
init_printing()
from IPython.display import display

# This command makes plots appear inside the browser window
%matplotlib inline
```

## Mass-spring-damper system

The differential equation that governs a single degree-of-freedom mass-spring-damper system, with a forcing term $f(t)$, is

$$
m \frac{d^{2}y}{dt^{2}} + \lambda \frac{dy}{dt} + ky = f(t)
$$

To solve this problem using SymPy, we first define the symbols $t$ (time), $m$ (mass), $\lambda$ (damper coefficient) and $k$ (spring stiffness), and the function $y$ (displacement): 


```
t, m, lmbda, k = symbols("t m lambda k")
y = Function("y")
f = Function("f")
```

Next, we define the differential equation, and print it to the screen:


```
eqn = Eq(m*Derivative(y(t), t, t) + lmbda*Derivative(y(t), t) + k*y(t), f(t))
display(eqn)
```

## Classifying the ODE

Checking the order of the ODE:


```
print("This order of the ODE is: {}".format(ode_order(eqn, y(t))))
```

    This order of the ODE is: 2


As expected, it is second order. Now classifying the ODE:


```
print("Properties of the ODE are: {}".format(classify_ode(eqn)))
```

    Properties of the ODE are: ('nth_linear_constant_coeff_variation_of_parameters', 'nth_linear_constant_coeff_variation_of_parameters_Integral')


## Solving the ODE

We can now try to solve the ODE symbolically, but without defining $f(t)$ the result will be umanageble. Sympy will return a solution but it will too complicated to interpret.

We will therefore specify some of the parameters in the equation:


```
m0, lmbda0, k0 = 4, 1, 16
eqn1 = eqn.subs(m, m0).subs(k, k0).subs(lmbda, lmbda0)

dsolve(eqn1)
```

The solution is looking more manageable, but we still haven't specified $f(t)$. Setting $f(t) = 2$ and then solving,


```
eqn1 = eqn1.subs(f(t), 2)
print("ODE with f(t) = 2")
display(eqn1)

ys = dsolve(eqn1, y(t))
print("Solution to ODE:")
display(ys)
```

The solution is now looking tractable. If we set $C_{1} = C_{2} = 1$, we can plot the solution. 


```
ys = ys.subs('C1', 1).subs('C2', 1)
plot(ys.args[1], (t, 0.0, 20.0), xlabel="time", ylabel="displacement")
```

We can see the that the solution is a decaying sinusoidal function.

### Exercise: 

Try changing the values for $m$, $\lambda$ and $k$ to investigate their influence in the response. You can try negative values, but only positive values are physically possible.

### Exercise:

Try changing $f(t)$ to explore the response for different forcing functions.

### Example: resonance

An undamped system ($\lambda = 0$) changes character when the forcing term is sinusoidal at the natural frequency of the system. The corresponds to the case of the forcing term having the same form as the solution to the homogeneous problem.

The natural frequency $\omega_{n} = \sqrt{k/m}$. Setting $f(t) = \sin( \omega x)$



```
# Define natural frequency
omega = Symbol("omega")
omega = sqrt(k/m)

# Set f(t) = sin(\omega t)
eqn1 = eqn.subs(f(t), sin(omega*t))
display(eqn1)
```

We now set values for $m$, $\lambda$ and $k$ (we'll set $\lambda = 0$ which corresponds to an undamped system):


```
# Set numerical values for m, lambda and k
m0, lmbda0, k0 = 4, 0, 1
eqn1 = eqn1.subs(m, m0).subs(k, k0).subs(lmbda, lmbda0)

# Solve ODE
display(eqn1)
y = dsolve(eqn1)
display(y)
```

Note the presence of $t$ in front of the cosine. This implies that the solution will grow linearly in time. We can see this by plotting the solution:


```
y = y.subs('C1', 1).subs('C2', 1)
plot(y.args[1], (t, 0.0, 100.0), xlabel="time", ylabel="displacement")
```

We observe that the solution grows in time - this is a feature of resonance.

### Example: Interactive control for the natural frequency of the mass-spring-damper system

####NOTE: This will not run on the online IPython Noteboook Viewer. You would need IPython 2.3.1 and run this cell locally for the interactive to work. 

The equation for the mass-spring-damper system above can be rewritten in terms of damping $\zeta$ and natural frequency $\omega_{n}$ as belowed:

$$ \frac{d^{2}y}{dt^{2}} + 2\zeta\omega_{n}\frac{dy}{dt} + \omega_{n}^2y = 0$$

which gives rise to three different solutions, whether the system is under-damped, critically damped or over-damped. 

The values of $\omega_{n}$ and $\zeta$ are defined as below, from the original mass-spring system:

$$ 
\omega_{n} = \sqrt{\frac{k}{m}}
$$
$$
\zeta = \frac{\lambda }{2 \sqrt{mk}}
$$

The interactive widget belows provides the visualisation between the changes of $m, k, \lambda$ on the final results.  



```
from sympy import *

# Import the interactive function from IPython widgets
from IPython.html.widgets import interact
from IPython.html import widgets 

# This initialises pretty printing
init_printing()
from IPython.display import display
%matplotlib inline 

# Declare modules nump
from numpy import *
from matplotlib.pyplot import *

# Definie time limits and step. This gives a time limit from 0 to 500 with a step of 5.
t = np.linspace(0,10,500)            

# Define the initial conditions velocity = 0, position = 1
y0 = [1.0, 0.0]

# Set up the function that plots the repsonse based on slider changes
# These values of m, k and lambda here are in SI units. m in kg, k in N/m and lambda in Ns/m

def plot_response(m = 5, k = 4, lmbda = 0.0):
   
    # Calculate and print out values of omega_n and zeta
    omega_n = sqrt(k/m)
    zeta = lmbda/(2*sqrt(m*k))
    print "Natural frequency is",omega_n
    print "Damping ratio is",zeta
    
    # Define natural and damped frequency
    omega_d = omega_n*sqrt(1-zeta**2)
    
    # Define y(t) based on different values of zeta
    
    if zeta < 1: # Under-damped
        A = y0[0];
        B = (1/omega_d)*(zeta*omega_n*y0[0] + y0[1])
        y = exp(-zeta*omega_n*t)*(A*cos(omega_d*t) + B*sin(omega_d*t))
    elif zeta == 1: # Critically damped
        A = y0[0]
        B = y0[1] + omega_n*y0[0]
        y = (A + B*t)*exp(-omega_n*t)
    elif zeta > 1: # Over-damped. root_p and root_n are the positive and negative solutions of the characteristic equations. 
        delta = (2*zeta*omega_n)**2 - 4*(omega_n**2)
        root_p = (0.5)*(-2*zeta*omega_n + sqrt(delta))
        root_n = (0.5)*(-2*zeta*omega_n - sqrt(delta))
        A = y0[0] + (root_p*y0[0] - y0[1])/(root_n - root_p)
        B = (-1)*(root_p*y0[0] - y0[1])/(root_n - root_p)
        y = A*exp(root_p*t) + B*exp(root_n*t)
        
    
    # Draw an label the figure window
    ax = gca()
    ax.grid(True,linestyle=':',color='0.5')
    xlabel('Time (s)',fontsize=20,labelpad=5)
    ylabel('Position (m)',fontsize=20,labelpad=10)
    ylim(-1.2,1.2)
    xlim(0,10)
    
    # Normalised value of y
    y = y/y[0]
    
    # Plot the function on the figure
    plot(t,y,color="red",linewidth=1.5)
   

# Call the slider interaction
interact(plot_response, m=(0,20, 0.1), k = (0,10,0.1), lmbda = (0,10,0.1))
```


```

```

2e7f38c68684f432a5e82f3d7f795d8c30dad26d
ipynb
Jupyter Notebook
notebooks/Lecture3.ipynb
quang-ha/IA-maths-Ipython
8ff8533d64a3d8db8e4813a7b6dfee39339fd846
BSD-3-Clause
notebooks/Lecture3.ipynb
quang-ha/IA-maths-Ipython
8ff8533d64a3d8db8e4813a7b6dfee39339fd846
BSD-3-Clause
notebooks/Lecture3.ipynb
quang-ha/IA-maths-Ipython
8ff8533d64a3d8db8e4813a7b6dfee39339fd846
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES