```python
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.integrate as inte
```

### ODE
$$F(x,y,y',...,y^{(n-1)})=y^{(n)}$$
where $y^{(n)}$ is the n-th derivative.
- Order of an ODE is the highest order derivtive that appears
- Mechanical systems are usually second order
- Recall Newton's second $\mathbf{F=ma=m\ddot{x}}$

### State-space for 2nd order systems
Define vector equation

$$x = \begin{pmatrix} x_1\\x_2\end{pmatrix}=\begin{pmatrix}\chi\\\dot{\chi}\end{pmatrix}\Longrightarrow \dot{x}=\begin{pmatrix}x_2\\\frac{F}{m}\end{pmatrix}$$

### Numerical ODE integration
- Consider $\dot{x}=\alpha$
- Could use fixed timestep, $dt$
- Set $x(t_{k+1})=x(t_k)+\alpha dt$
- If $\alpha(t)$ is not fixed, the results will be inaccurate
- ode estimates how the right side is changing and picks the best timestep $dt$.
- Only works when right hand side is smooth

### Examples from MATLAB to python
Try to solve the following ode problems from https://uk.mathworks.com/help/matlab/ref/ode45.html using python.

### Solve ODE with Single Solution Component
$$y'=2t$$


```python
def odeSSC(t,y):
    return 2*t

t0, t1 = 0, 5
t = np.linspace(t0,t1,100)
y0 = 0
y = np.zeros((len(t),1))
y[0,:] = y0 
r = inte.ode(odeSSC).set_integrator("dopri5")  # choice of method
r.set_initial_value(y0, t0)   # initial values
for i in range(1, t.size):
   y[i, :] = r.integrate(t[i]) # get one more value, add it to the array
   if not r.successful():
       raise RuntimeError("Could not integrate")
plt.plot(t, y)
```

### Solve Nonstiff equation
The van der Pol equation is a second order ODE:
$$y_1''-\mu(1-y_1^2)y_1'+y_1=0$$
where $\mu>0$ is a scalar parameter. Rewrite this equation as first order ODEs by making the substitution $y_2=y_1'$. The resulting system of first-order ODEs is, 

\begin{align}
y_1' &= y_2 \\
y_2' &= \mu(1-y_1^2)y_2-y_1
\end{align}


```python
def vdp1(t, y):
    return np.array([y[1], (1 - y[0]**2)*y[1] - y[0]])

t0, t1 = 0, 20                # start and end
t = np.linspace(t0, t1, 100)  # the points of evaluation of solution
y0 = [2, 0]                   # initial value
y = np.zeros((len(t), len(y0)))   # array for solution
y[0, :] = y0
r = inte.ode(vdp1).set_integrator("dopri5")  # choice of method
r.set_initial_value(y0, t0)   # initial values
for i in range(1, t.size):
   y[i, :] = r.integrate(t[i]) # get one more value, add it to the array
   if not r.successful():
       raise RuntimeError("Could not integrate")
plt.plot(t, y)
```

### Pass Extra Parameters to ODE Function
Solve the ODE,
$$y''=\frac{A}{B}ty$$
Rewriting the equation as a 1st order ODE yields,

\begin{align}
y_1' &= y_2\\
y_2' &= \frac{A}{B}ty_1
\end{align}


```python
def odefcn(t,y,A,B):
    return np.array([y[1],A/B*t*y[0]])
A,B = 1,2
t0, t1 = 0, 5
t = np.linspace(t0,t1,100)
y0 = [0,0.01]                   # initial value
y = np.zeros((len(t), len(y0)))   # array for solution
r = inte.ode(lambda t,y:odefcn(t,y,A,B)).set_integrator("dopri5")  # choice of method
r.set_initial_value(y0, t0)   # initial values
for i in range(1, t.size):
   y[i, :] = r.integrate(t[i]) # get one more value, add it to the array
   if not r.successful():
       raise RuntimeError("Could not integrate")
plt.plot(t, y)
```

### ODE with Time-Dependent Terms
Consider the following ODE with time-dependent parameters,
$$y'(t)+f(t)y(t) = g(t)$$


```python
ft = np.linspace(0,5,25)
gt = np.linspace(1,6,25)
f = ft**2-ft-3
g = 3*np.sin(gt-.25)
t0, t1, ic, y0 = 1, 5, 1, 1
t= np.linspace(t0,t1,100)
```


```python
def myode(t,y,ft,f,gt,g):
    f = np.interp(t,ft,f)
    g = np.interp(t,gt,g)
    dydt = -f*y+g
    return dydt
```


```python
y = np.zeros((len(t), 1))   # array for solution
r = inte.ode(lambda t,y:myode(t,y,ft,f,gt,g)).set_integrator("dopri5",atol=1e-4,rtol=1e-2)  # choice of method
r.set_initial_value(y0, t0)   # initial values
for i in range(1, t.size):
    tmp = r.integrate(t[i]) # get one more value, add it to the array
    y[i, :] = tmp
    if not r.successful():
        raise RuntimeError("Could not integrate")
plt.plot(t, y)
```

### Evaluate and Extend Solution Structure
The van der Pol equation is a second order ODE
$$y_1''-\mu(1-y_1^2)y_1'+y_1=0$$


```python
t0, t1 = 0, 20                # start and end
t = np.linspace(t0, t1, 100)  # the points of evaluation of solution
y0 = [2, 0]                   # initial value
y = np.zeros((len(t), len(y0)))   # array for solution
y[0, :] = y0
r = inte.ode(vdp1).set_integrator("dopri5")  # choice of method
r.set_initial_value(y0, t0)   # initial values
for i in range(1, t.size):
   y[i, :] = r.integrate(t[i]) # get one more value, add it to the array
   if not r.successful():
       raise RuntimeError("Could not integrate")
```


```python
plt.plot(t,y[:,0])
```


```python
x0, x1 = 20,35
x = np.linspace(x0,x1,350)
yx0 = y[-1,:]
yx = np.zeros((len(x), len(yx0)))   # array for solution
yx[0,:] = yx0

r.set_initial_value(yx0, x0)   # initial values
for i in range(1, x.size):
   yx[i, :] = r.integrate(x[i]) # get one more value, add it to the array
   if not r.successful():
       raise RuntimeError("Could not integrate")
    
plt.plot(t,y[:,0])
plt.plot(x,yx[:,0])
```

### Robotics Capstone 
#### 1. Solving ODE
$$\dot{x}=\begin{pmatrix}\dot{x_1}\\\dot{x_2}\end{pmatrix}=\begin{pmatrix} x_2 \\ \gamma\cos(\omega t)-\alpha x_1-\delta x_2 -\beta x_1^3\end{pmatrix}$$


```python
def dyn(t,X,beta):
    delta,alpha,gamma,omega = 1,1,1,1
    x = X[0]; xd = X[1]
    return np.array([xd,gamma*np.cos(omega*t)-alpha*x-delta*xd-beta*x**3])
```


```python
def simode(beta,tend):
    t0, t1 = 0, tend                # start and end
    t = np.linspace(t0, t1, 100)  # the points of evaluation of solution
    y0 = [1, 1]                   # initial value
    y = np.zeros((len(t), len(y0)))   # array for solution
    y[0, :] = y0
    r = inte.ode(lambda t,X: dyn(t,X,beta)).set_integrator("dopri5")  
    r.set_initial_value(y0, t0)   # initial values
    for i in range(1, t.size):
       y[i, :] = r.integrate(t[i]) # get one more value, add it to the array
       if not r.successful():
           raise RuntimeError("Could not integrate")
    
    return t,y
```


```python
t,y = simode(1,50)
plt.figure(figsize=(20,5))
plt.subplot(1,2,1); plt.plot(t,y)
plt.subplot(1,2,2); plt.plot(y[:,0], y[:,1])
```


```python
t,y = simode(5,50)
plt.figure(figsize=(20,5))
plt.subplot(1,2,1); plt.plot(t,y)
plt.subplot(1,2,2); plt.plot(y[:,0], y[:,1])
```


```python
t,y = simode(10,50)
plt.figure(figsize=(20,5))
plt.subplot(1,2,1); plt.plot(t,y)
plt.subplot(1,2,2); plt.plot(y[:,0], y[:,1])
```


```python

```

6b84ddfebcb372591419b45fa23a0f74d1869f1f
ipynb
Jupyter Notebook
Robotics/ODEsolver/.ipynb_checkpoints/Ordinary Differential Equation-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2019-11-20T10:20:29.000Z
2021-11-09T11:15:23.000Z
Capstone/ODE solver/Ordinary Differential Equation.ipynb
kasiv008/Robotics
302b3336005acd81202ebbbb0c52a4b2692fa9c7
MIT
Capstone/ODE solver/Ordinary Differential Equation.ipynb
kasiv008/Robotics
302b3336005acd81202ebbbb0c52a4b2692fa9c7
MIT
2020-05-27T03:56:38.000Z
2021-05-02T13:15:42.000Z
Qwen/Qwen-72B
1. YES
2. YES