# Solving ODEs with scipy.integrate.solve_ivp

## Solving ordinary differential equations (ODEs)

Here we will revisit the differential equations solved in 5300_Jupyter_Python_intro_01.ipynb with `odeint`, only now we'll use `solve_ivp` from Scipy.  We'll compare the new and old solutions as we go.

### First-order ODE


```python
# Import the required modules
import numpy as np
import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp   # Now preferred to odeint
```

Let's try a one-dimensional first-order ODE, say:

$\begin{align}
\quad 
\frac{dv}{dt} = -g, \quad \mbox{with} \quad v(0) = 10
\end{align}$

in some appropriate units (we'll use MKS units by default).  This ODE can be separated and directly integrated:

$\begin{align}
  \int_{v_0=10}^{v} dv' = - g \int_{0}^{t} dt'
  \quad\Longrightarrow\quad
    v - v_0 = - g (t - 0)
  \quad\Longrightarrow\quad
   v(t) = 10 - gt
\end{align}$



The goal is to find the solution $v(t)$ as an array `v_pts` at the times in the array `t_pts`.


```python
# Define a function which calculates the derivative
def dv_dt_new(t, v, g=9.8):
    """Returns the right side of a simple first-order ODE with default g."""
    return -g   

t_start = 0.
t_end = 10.
t_pts = np.linspace(t_start, t_end, 20)  # 20 points between t=0 and t=10.

v_0 = np.array([10.0])  # initial condition, in form of a list or numpy array

abserr = 1.e-8
relerr = 1.e-8

solution = solve_ivp(dv_dt_new, (t_start, t_end), v_0, t_eval=t_pts,
                     rtol=relerr, atol=abserr)  
    # solve_ivp( function for rhs with (t, v) argument (cf. (v,t) for odeint), 
    #            tspan=(starting t value, ending t value),
    #            initial value of v(t), array of points we want to know v(t),
    #            method='RK45' is the default method,
    #            rtol=1.e-3, atol=1.e-6 are default tolerances
    #          )
v_pts = solution.y  # array of results at t_pts
```


```python
v_pts.shape   # 1 x 100 matrix (row vector)
```




    (1, 20)



Here's how we did it before with odeint:


```python
from scipy.integrate import odeint   

# Define a function which calculates the derivative
def dv_dt(v, t, g=9.8):
    """Returns the right side of a simple first-order ODE with default g."""
    return -g   

t_pts = np.linspace(0., 10., 20)     # 20 points between t=0 and t=10.
v_0 = 10.  # the initial condition
v_pts_odeint = odeint(dv_dt, v_0, t_pts)  # odeint( function for rhs, 
                                          #         initial value of v(t),
                                          #         array of t values )

```


```python
v_pts_odeint.shape   # 100 x 1 matrix (column vector)
```




    (20, 1)



Make a table comparing results (using `flatten()` to make the matrices into arrays):


```python
print('    t     v(t) [solve_ivp]    v(t) [odeint]')
for t, v_solve_ivp, v_odeint in zip(t_pts, 
                                    v_pts.flatten(), 
                                    v_pts_odeint.flatten()):
    print(f' {t:6.3f}   {v_solve_ivp:12.7f}       {v_odeint:12.7f}')
```

        t     v(t) [solve_ivp]    v(t) [odeint]
      0.000     10.0000000         10.0000000
      0.526      4.8421053          4.8421053
      1.053     -0.3157895         -0.3157895
      1.579     -5.4736842         -5.4736842
      2.105    -10.6315789        -10.6315789
      2.632    -15.7894737        -15.7894737
      3.158    -20.9473684        -20.9473684
      3.684    -26.1052632        -26.1052632
      4.211    -31.2631579        -31.2631579
      4.737    -36.4210526        -36.4210526
      5.263    -41.5789474        -41.5789474
      5.789    -46.7368421        -46.7368421
      6.316    -51.8947368        -51.8947368
      6.842    -57.0526316        -57.0526316
      7.368    -62.2105263        -62.2105263
      7.895    -67.3684211        -67.3684211
      8.421    -72.5263158        -72.5263158
      8.947    -77.6842105        -77.6842105
      9.474    -82.8421053        -82.8421053
     10.000    -88.0000000        -88.0000000


Differences between `solve_ivp` and `odeint`:
* `dv_dt(t, v)`  vs.  `dv_dt(v, t)`, i.e., the function definitions have the arguments reversed.
* With `odeint`, you only specify the full array of $t$ points you want to know $v(t)$ at.  With `solve_ivp`, you first specify the starting $t$ and ending $t$ as a tuple: `(t_start, t_end)` and then (optionally) specify `t_eval=t_pts` to evaluate $v$ at the points in the `t_pts` array.
* `solve_ivp` returns an object from which $v(t)$ (and other results) can be found, while `ode_int` returns $v(t)$.
* For this single first-order equation, $v(t)$ is returned for the $N$ requested $t$ points as a $1 \times N$ two-dimensional array by `solve_ivp` and as a $N \times 1$ array by `odeint`.
* `odeint` has no choice of solver while the `solve_ivp` solver can be set by `method`.  The default is `method='RK45'`, which is good, general-purpose Runge-Kutta solver.  

### Second-order ODE

Suppose we have a second-order ODE such as:

$$
\quad y'' + 2 y' + 2 y = \cos(2x), \quad \quad y(0) = 0, \; y'(0) = 0
$$

We can turn this into two first-order equations by defining a new dependent variable. For example,

$$
\quad z \equiv y' \quad \Rightarrow \quad z' + 2 z + 2y = \cos(2x), \quad z(0)=y(0) = 0.
$$

Now introduce the vector 

$$
  \mathbf{U}(x) = \left(\begin{array}{c}
                         y(x) \\
                         z(x)
                        \end{array}
                  \right)
        \quad\Longrightarrow\quad
    \frac{d\mathbf{U}}{dx} = \left(\begin{array}{c}
                                    z \\
                                    -2 y' - 2 y + \cos(2x)
                                   \end{array}
                             \right) 
$$

We can solve this system of ODEs using `solve_ivp` with lists, as follows.  We will try it first without specifying the relative and absolute error tolerances rtol and atol.


```python
# Define a function for the right side
def dU_dx_new(x, U):
    """Right side of the differential equation to be solved.
    U is a two-component vector with y=U[0] and z=U[1]. 
    Thus this function should return [y', z']
    """
    return [U[1], -2*U[1] - 2*U[0] + np.cos(2*x)]

# initial condition U_0 = [y(0)=0, z(0)=y'(0)=0]
U_0 = [0., 0.]

x_pts = np.linspace(0, 15, 20)  # Set up the mesh of x points
result = solve_ivp(dU_dx_new, (0, 15), U_0, t_eval=x_pts)
y_pts = result.y[0,:]   # Ok, this is tricky.  For each x, result.y has two 
                        #  components.  We want the first component for all
                        #  x, which is y(x).  The 0 means the first index and 
                        #  the : means all of the x values.

```

Here's how we did it before with `odeint`:


```python
# Define a function for the right side
def dU_dx(U, x):
    """Right side of the differential equation to be solved.
    U is a two-component vector with y=U[0] and z=U[1]. 
    Thus this function should return [y', z']
    """
    return [U[1], -2*U[1] - 2*U[0] + np.cos(2*x)]

# initial condition U_0 = [y(0)=0, z(0)=y'(0)=0]
U_0 = [0., 0.]

x_pts = np.linspace(0, 15, 20)  # Set up the mesh of x points
U_pts = odeint(dU_dx, U_0, x_pts)  # U_pts is a 2-dimensional array
y_pts_odeint = U_pts[:,0]  # Ok, this is tricky.  For each x, U_pts has two 
                           #  components.  We want the upper component for all
                           #  x, which is y(x).  The : means all of the first 
                           #  index, which is x, and the 0 means the first
                           #  component in the other dimension.
```

Make a table comparing results (using `flatten()` to make the matrices into arrays):


```python
print('    x     y(x) [solve_ivp]    y(x) [odeint]')
for x, y_solve_ivp, y_odeint in zip(x_pts, 
                                    y_pts.flatten(), 
                                    y_pts_odeint.flatten()):
    print(f' {x:6.3f}   {y_solve_ivp:12.7f}       {y_odeint:12.7f}')
```

        x     y(x) [solve_ivp]    y(x) [odeint]
      0.000      0.0000000          0.0000000
      0.789      0.1360331          0.1360684
      1.579      0.0346990          0.0347028
      2.368     -0.2285865         -0.2287035
      3.158     -0.0975122         -0.0974702
      3.947      0.2065834          0.2067492
      4.737      0.0927010          0.0927536
      5.526     -0.2041225         -0.2042677
      6.316     -0.0865181         -0.0865921
      7.105      0.2064689          0.2066669
      7.895      0.0832986          0.0832707
      8.684     -0.2080277         -0.2081975
      9.474     -0.0799116         -0.0799972
     10.263      0.2092969          0.2094602
     11.053      0.0765611          0.0765810
     11.842     -0.2104970         -0.2107011
     12.632     -0.0731547         -0.0731411
     13.421      0.2117489          0.2118952
     14.211      0.0695941          0.0696868
     15.000     -0.2129369         -0.2130316


Not very close agreement by the end.  Run both again with greater accuracy.


```python
relerr = 1.e-10
abserr = 1.e-10

result = solve_ivp(dU_dx_new, (0, 15), U_0, t_eval=x_pts, 
                   rtol=relerr, atol=abserr)
y_pts = result.y[0,:]    

U_pts = odeint(dU_dx, U_0, x_pts, 
               rtol=relerr, atol=abserr)  
y_pts_odeint = U_pts[:,0]   

print('    x     y(x) [solve_ivp]    y(x) [odeint]')
for x, y_solve_ivp, y_odeint in zip(x_pts, 
                                    y_pts.flatten(), 
                                    y_pts_odeint.flatten()):
    print(f' {x:6.3f}   {y_solve_ivp:12.7f}       {y_odeint:12.7f}')
```

        x     y(x) [solve_ivp]    y(x) [odeint]
      0.000      0.0000000          0.0000000
      0.789      0.1360684          0.1360684
      1.579      0.0347028          0.0347028
      2.368     -0.2287035         -0.2287035
      3.158     -0.0974702         -0.0974702
      3.947      0.2067492          0.2067492
      4.737      0.0927536          0.0927536
      5.526     -0.2042678         -0.2042678
      6.316     -0.0865921         -0.0865921
      7.105      0.2066669          0.2066669
      7.895      0.0832707          0.0832707
      8.684     -0.2081975         -0.2081975
      9.474     -0.0799972         -0.0799972
     10.263      0.2094602          0.2094602
     11.053      0.0765810          0.0765810
     11.842     -0.2107011         -0.2107011
     12.632     -0.0731411         -0.0731411
     13.421      0.2118952          0.2118952
     14.211      0.0696868          0.0696868
     15.000     -0.2130316         -0.2130316


Comparing the results from when we didn't specify the errors we see that the default error tolerances for solve_ivp were insufficient.  Moral: specify them explicitly.  

a1167fae981759fdd686bf60c0bad38d601c2559
ipynb
Jupyter Notebook
week_4/ODEs_with_solve_ivp.ipynb
CLima86/Physics_5300_CDL
d9e8ee0861d408a85b4be3adfc97e98afb4a1149
MIT
week_4/ODEs_with_solve_ivp.ipynb
CLima86/Physics_5300_CDL
d9e8ee0861d408a85b4be3adfc97e98afb4a1149
MIT
week_4/ODEs_with_solve_ivp.ipynb
CLima86/Physics_5300_CDL
d9e8ee0861d408a85b4be3adfc97e98afb4a1149
MIT
Qwen/Qwen-72B
1. YES
2. YES