```python
%matplotlib inline
from matplotlib import pyplot as plt

import sympy as sp
import numpy as np
import math
```

### Exact solution using Sympy


```python
x = sp.Function("x")
K = sp.Function("K")
t = sp.symbols("t")
m = sp.symbols("m")
r = sp.symbols("r")
k = sp.symbols("k")

dgl = sp.Eq(m*x(t).diff(t).diff(t) + r*x(t).diff(t) + k*x(t), K(t))
dgl
```




$\displaystyle k x{\left(t \right)} + m \frac{d^{2}}{d t^{2}} x{\left(t \right)} + r \frac{d}{d t} x{\left(t \right)} = K{\left(t \right)}$




```python
sp.dsolve(dgl, x(t))
```




$\displaystyle x{\left(t \right)} = C_{1} e^{\frac{t \left(- r - \sqrt{- 4 k m + r^{2}}\right)}{2 m}} + C_{2} e^{\frac{t \left(- r + \sqrt{- 4 k m + r^{2}}\right)}{2 m}} - \frac{e^{\frac{t \left(- r - \sqrt{- 4 k m + r^{2}}\right)}{2 m}} \int K{\left(t \right)} e^{\frac{r t}{2 m}} e^{\frac{t \sqrt{- 4 k m + r^{2}}}{2 m}}\, dt}{\sqrt{- 4 k m + r^{2}}} + \frac{e^{\frac{t \left(- r + \sqrt{- 4 k m + r^{2}}\right)}{2 m}} \int K{\left(t \right)} e^{\frac{r t}{2 m}} e^{- \frac{t \sqrt{- 4 k m + r^{2}}}{2 m}}\, dt}{\sqrt{- 4 k m + r^{2}}}$



### Runge-Kutta solution using Scipy


```python
def oscillator(t, x, K_0, omega, m, r, k):

    x_1, x_2 = x # x: vector of x_1 and x_2
    
    dxdt = [x_2, 1/m*(K_0*math.sin(omega*t) - r*x_2 - k*x_1)]
    return dxdt
```


```python
from scipy.integrate import solve_ivp
init = [0, 0]
time = np.arange(0, 3, 0.1)
time_span = [0, 3]

K_0 = 1
omega = 0.2
m = 0.5
r = 0.1
k = 0.2

dgl_RK45 = solve_ivp(oscillator, time_span, init, method="RK45", t_eval=time, args = (K_0, omega, m, r, k))
```


```python
plt.figure()
plt.plot(time, dgl_RK45.y[0, :], dgl_RK45.y[1, :])
plt.legend(["$x_1$ Runge-Kutta", "$x_2$ Runge-Kutta"])
plt.xlabel("Time")
plt.ylabel("$y$")
plt.show()
```


```python
dgl_RK45
```




      message: 'The solver successfully reached the end of the integration interval.'
         nfev: 50
         njev: 0
          nlu: 0
          sol: None
       status: 0
      success: True
            t: array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2,
           1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2, 2.3, 2.4, 2.5,
           2.6, 2.7, 2.8, 2.9])
     t_events: None
            y: array([[0.00000000e+00, 6.63200869e-05, 5.27309831e-04, 1.76946532e-03,
            4.16797423e-03, 8.08480963e-03, 1.38672118e-02, 2.18506552e-02,
            3.23538823e-02, 4.56764372e-02, 6.20986655e-02, 8.18806121e-02,
            1.05259336e-01, 1.32462317e-01, 1.63694785e-01, 1.99137417e-01,
            2.38946337e-01, 2.83253116e-01, 3.32164771e-01, 3.85763766e-01,
            4.44108756e-01, 5.07242177e-01, 5.75185763e-01, 6.47938342e-01,
            7.25476244e-01, 8.07753300e-01, 8.94700842e-01, 9.86227707e-01,
            1.08222023e+00, 1.18254226e+00],
           [0.00000000e+00, 1.98600542e-03, 7.88292551e-03, 1.75874367e-02,
            3.09817736e-02, 4.79344570e-02, 6.83013237e-02, 9.19245793e-02,
            1.18635561e-01, 1.48256051e-01, 1.80598281e-01, 2.15465907e-01,
            2.52657766e-01, 2.91958231e-01, 3.33147574e-01, 3.76003879e-01,
            4.20303040e-01, 4.65818770e-01, 5.12322591e-01, 5.59583840e-01,
            6.07371908e-01, 6.55464606e-01, 7.03627220e-01, 7.51628873e-01,
            7.99245629e-01, 8.46260496e-01, 8.92463425e-01, 9.37651308e-01,
            9.81627978e-01, 1.02420421e+00]])
     y_events: None



145a9622960d3a86d6acfc335f73e5d685f50c2d
ipynb
Jupyter Notebook
04_mini-projects/00_mechanical_oscillator.ipynb
benebrue/py-algorithms-4-automotive-engineering
eb4c137f3125e3f758bd25bffc4f2e9e5800b517
MIT
04_mini-projects/00_mechanical_oscillator.ipynb
benebrue/py-algorithms-4-automotive-engineering
eb4c137f3125e3f758bd25bffc4f2e9e5800b517
MIT
04_mini-projects/00_mechanical_oscillator.ipynb
benebrue/py-algorithms-4-automotive-engineering
eb4c137f3125e3f758bd25bffc4f2e9e5800b517
MIT
Qwen/Qwen-72B
1. YES
2. YES