# ODEs (Lecture 7)


```python
# This cell just imports the relevant modules

import numpy as np
from math import pi, exp
from sympy import init_printing, sin, cos, Function, Symbol, diff, integrate, dsolve, checkodesol, solve, ode_order, classify_ode, pprint
import mpmath
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML
```

## Order of an ODE

**Slide 9**

Use `sympy` to define dependent and independent variables, constants, ODE, and to find the order of ODEs.


```python
t = Symbol('t')  # Independent variable
eta = Symbol('eta')  # Constant
v = Function('v')(t)  # Dependent variable v(t)
ode = diff(v,t) + eta*v # The ODE we wish to solve. Make sure the RHS is equal to zero.
print("ODE #1:") 
pprint(ode)
print("The order of ODE #1 is", ode_order(ode, v)) 

x = Function('x')(t)  # Dependent variable x(t)
m = Symbol('m')  # Constant
k = Symbol('k')  # Constant
ode = m*diff(x,t,2) + k*x
print("ODE #2:") 
pprint(ode)
print("The order of ODE #2 is", ode_order(ode, x)) 

y = Function('y')(t)  # Dependent variable y(t)
ode = diff(y,t,4) - diff(y,t,2)
print("ODE #3:") 
pprint(ode)
print("The order of ODE #3 is", ode_order(ode, y))
```

    ODE #1:
             d       
    η⋅v(t) + ──(v(t))
             dt      
    The order of ODE #1 is 1
    ODE #2:
                 2      
                d       
    k⋅x(t) + m⋅───(x(t))
                 2      
               dt       
    The order of ODE #2 is 2
    ODE #3:
        2           4      
       d           d       
    - ───(y(t)) + ───(y(t))
        2           4      
      dt          dt       
    The order of ODE #3 is 4


# Analytical solutions

**Slide 14**

Solving ODEs analytically using `sympy.dsolve`


```python
x = Symbol('x')  # Independent variable
y = Function('y')(x)  # Dependent variable y(x)

# The ODE we wish to solve. Make sure the RHS is equal to zero.
ode = diff(y,x) - 2*x*(1-y)
solution = dsolve(ode, y) # Solve the ode for function y(x).
print("ODE #4:") 
pprint(ode)
print("The solution to ODE #4 is: ", solution) 
```

    ODE #4:
                      d       
    -2⋅x⋅(1 - y(x)) + ──(y(x))
                      dx      
    The solution to ODE #4 is:  Eq(y(x), C1*exp(-x**2) + 1)



```python
x_3d = np.arange(-5, 5, 0.01)
y_3d = np.arange(-5, 5, 0.01)
X, Y = np.meshgrid(x_3d, y_3d)

dydx = 2 * X * (1-Y)

x = np.linspace(-5, 5, 1000)  
y = 1 + np.exp(-x**2)

fig = plt.figure(figsize=(16,8))

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(X, Y, dydx, cmap='seismic', edgecolor='k', lw=0.25)
ax1.set_xlabel('x')
ax1.set_ylabel('y(x)')
ax1.set_zlabel('dy/dx')
ax1.set_title('dy/dx = 2x(1-y)', fontsize=14)

ax2 = fig.add_subplot(122)
ax2.plot(x, y, 'k', label='y(x) = 1 + exp(-x**2)')
ax2.set_xlabel('x')
ax2.set_ylabel('y(x)')
ax2.set_title("Solution to ODE 4", fontsize=14)
ax2.legend(loc='best')
ax2.grid(True)

fig.tight_layout()
plt.show()
```

```{note}

The function `checkodesol` checks that the result from `dsolve` is indeed a solution to the ode. It substitutes in 'solution' into 'ode' and checks that the RHS is zero. If it is, the function returns 'True'.

```


```python
print("Checking solution using checkodesol...") 
check = checkodesol(ode, solution)

print("Output from checkodesol:", check)

if(check[0] == True):
    print("y(x) is indeed a solution to ODE #4") 
else:
    print("y(x) is NOT a solution to ODE #4") 
```

    Checking solution using checkodesol...
    Output from checkodesol: (True, 0)
    y(x) is indeed a solution to ODE #4


```{note}

The mpmath module can handle initial conditions (x0, y0) when solving an initial value problem, using the 
odefun function. However, this will not give you an analytical solution to the ODE, only a numerical solution.

```


```python
f = mpmath.odefun(lambda x, y: 2*x*(1-y), x0=0, y0=2)

# compares the numerical solution f(x) with the values of the (already known) analytical solution 
# between x=0 and x=10
for x in np.linspace(0, 10, 101):
    print("x=%.1f" % (x), ",", f(x), ",", 1+exp(-x**2))
```

    x=0.0 , 2.0 , 2.0
    x=0.1 , 1.99004983374917 , 1.990049833749168
    x=0.2 , 1.96078943915232 , 1.960789439152323
    x=0.3 , 1.91393118527123 , 1.9139311852712282
    x=0.4 , 1.85214378896621 , 1.8521437889662113
    x=0.5 , 1.7788007830714 , 1.778800783071405
    x=0.6 , 1.69767632607103 , 1.697676326071031
    x=0.7 , 1.61262639418442 , 1.612626394184416
    x=0.8 , 1.52729242404305 , 1.5272924240430485
    x=0.9 , 1.44485806622294 , 1.444858066222941
    x=1.0 , 1.36787944117144 , 1.3678794411714423
    x=1.1 , 1.29819727942989 , 1.2981972794298873
    x=1.2 , 1.23692775868212 , 1.2369277586821217
    x=1.3 , 1.18451952399299 , 1.1845195239929893
    x=1.4 , 1.14085842092104 , 1.1408584209210448
    x=1.5 , 1.10539922456186 , 1.1053992245618642
    x=1.6 , 1.0773047404433 , 1.0773047404432998
    x=1.7 , 1.05557621261148 , 1.055576212611483
    x=1.8 , 1.03916389509899 , 1.039163895098987
    x=1.9 , 1.02705184686635 , 1.0270518468663503
    x=2.0 , 1.01831563888873 , 1.0183156388887342
    x=2.1 , 1.01215517832991 , 1.012155178329915
    x=2.2 , 1.00790705405159 , 1.0079070540515935
    x=2.3 , 1.00504176025969 , 1.005041760259691
    x=2.4 , 1.00315111159844 , 1.0031511115984444
    x=2.5 , 1.00193045413623 , 1.0019304541362277
    x=2.6 , 1.0011592291739 , 1.0011592291739047
    x=2.7 , 1.00068232805276 , 1.0006823280527564
    x=2.8 , 1.00039366904066 , 1.0003936690406552
    x=2.9 , 1.00022262985692 , 1.0002226298569188
    x=3.0 , 1.00012340980409 , 1.0001234098040868
    x=3.1 , 1.0000670548243 , 1.000067054824303
    x=3.2 , 1.00003571284964 , 1.0000357128496415
    x=3.3 , 1.00001864374233 , 1.0000186437423315
    x=3.4 , 1.00000954016287 , 1.000009540162873
    x=3.5 , 1.00000478511739 , 1.000004785117392
    x=3.6 , 1.0000023525752 , 1.0000023525752
    x=3.7 , 1.00000113372714 , 1.0000011337271388
    x=3.8 , 1.00000053553478 , 1.0000005355347803
    x=3.9 , 1.0000002479596 , 1.0000002479596017
    x=4.0 , 1.00000011253517 , 1.0000001125351747
    x=4.1 , 1.00000005006218 , 1.0000000500621802
    x=4.2 , 1.00000002182958 , 1.000000021829578
    x=4.3 , 1.00000000933029 , 1.0000000093302877
    x=4.4 , 1.00000000390894 , 1.0000000039089385
    x=4.5 , 1.00000000160523 , 1.000000001605228
    x=4.6 , 1.00000000064614 , 1.0000000006461431
    x=4.7 , 1.00000000025494 , 1.0000000002549383
    x=4.8 , 1.0000000000986 , 1.0000000000985951
    x=4.9 , 1.00000000003738 , 1.0000000000373757
    x=5.0 , 1.00000000001389 , 1.000000000013888
    x=5.1 , 1.00000000000506 , 1.0000000000050582
    x=5.2 , 1.00000000000181 , 1.0000000000018059
    x=5.3 , 1.00000000000063 , 1.000000000000632
    x=5.4 , 1.00000000000022 , 1.0000000000002167
    x=5.5 , 1.00000000000007 , 1.0000000000000728
    x=5.6 , 1.00000000000002 , 1.000000000000024
    x=5.7 , 1.00000000000001 , 1.0000000000000078
    x=5.8 , 1.0 , 1.0000000000000024
    x=5.9 , 1.0 , 1.0000000000000007
    x=6.0 , 1.0 , 1.0000000000000002
    x=6.1 , 1.0 , 1.0
    x=6.2 , 1.0 , 1.0
    x=6.3 , 1.0 , 1.0
    x=6.4 , 1.0 , 1.0
    x=6.5 , 1.0 , 1.0
    x=6.6 , 1.0 , 1.0
    x=6.7 , 1.0 , 1.0
    x=6.8 , 1.0 , 1.0
    x=6.9 , 1.0 , 1.0
    x=7.0 , 1.0 , 1.0
    x=7.1 , 1.0 , 1.0
    x=7.2 , 1.0 , 1.0
    x=7.3 , 1.0 , 1.0
    x=7.4 , 1.0 , 1.0
    x=7.5 , 1.0 , 1.0
    x=7.6 , 1.0 , 1.0
    x=7.7 , 1.0 , 1.0
    x=7.8 , 1.0 , 1.0
    x=7.9 , 1.0 , 1.0
    x=8.0 , 1.0 , 1.0
    x=8.1 , 1.0 , 1.0
    x=8.2 , 1.0 , 1.0
    x=8.3 , 1.0 , 1.0
    x=8.4 , 1.0 , 1.0
    x=8.5 , 1.0 , 1.0
    x=8.6 , 1.0 , 1.0
    x=8.7 , 1.0 , 1.0
    x=8.8 , 1.0 , 1.0
    x=8.9 , 1.0 , 1.0
    x=9.0 , 1.0 , 1.0
    x=9.1 , 1.0 , 1.0
    x=9.2 , 1.0 , 1.0
    x=9.3 , 1.0 , 1.0
    x=9.4 , 1.0 , 1.0
    x=9.5 , 1.0 , 1.0
    x=9.6 , 1.0 , 1.0
    x=9.7 , 1.0 , 1.0
    x=9.8 , 1.0 , 1.0
    x=9.9 , 1.0 , 1.0
    x=10.0 , 1.0 , 1.0


# Separation of variables

**Slide 20**

We can solve ODEs via separation of variables in Python using `sympy.dsolve` by passing the `hint` argument. Note that the optional `hint` argument here has been used to tell SymPy how to solve the ODE. However, it is usually smart enough to work it out for itself.


```python
x = Symbol('x')  # Independent variable
y = Function('y')(x)  # Dependent variable y(x)

# The ODE we wish to solve.
ode = (1.0/y)*diff(y,x) - cos(x)
print("ODE #5:") 
pprint(ode)

# Solve the ode for function y(x).using separation of variables.
solution = dsolve(ode, y, hint='separable')
print("The solution to ODE #5 is:") 
pprint(solution)
```

    ODE #5:
                  d       
              1.0⋅──(y(x))
                  dx      
    -cos(x) + ────────────
                  y(x)    
    The solution to ODE #5 is:
               1.0⋅sin(x)
    y(x) = C₁⋅ℯ          



```python
x_3d = np.arange(-5, 5, 0.01)
y_3d = np.arange(-5, 5, 0.01)
X, Y = np.meshgrid(x_3d, y_3d)

dydx = Y * np.cos(X) 

x = np.linspace(-5, 5, 1000)
y = np.exp(np.sin(x))
    
fig = plt.figure(figsize=(16, 8))

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(X, Y, dydx, cmap='seismic', edgecolor='k', lw=0.25)
ax1.set_xlabel('x')
ax1.set_ylabel('y(x)')
ax1.set_zlabel('dy/dx')
ax1.set_title('dy/dx = ycos(x)', fontsize=14)

ax2 = fig.add_subplot(122)
ax2.plot(x, y, 'k', label='y(x) = exp(sin(x))')
ax2.set_xlabel('x')
ax2.set_ylabel('y(x)')
ax2.set_title("Solution of ODE 5", fontsize=14)
ax2.legend(loc='best')
ax2.grid(True)

fig.tight_layout()
plt.show()
```

# Integration factor

**Slide 23**


```python
x = Symbol('x')  # Independent variable
y = Function('y')(x)  # Dependent variable y(x)
# The ODE we wish to solve.
ode = diff(y,x) - 2*x + 2*x*y

print("ODE #6:") 
pprint(ode)
# Solve the ode for function y(x).using separation of variables
solution = dsolve(ode, y)
print("The solution to ODE #6 is:", solution) 
```

    ODE #6:
                     d       
    2⋅x⋅y(x) - 2⋅x + ──(y(x))
                     dx      
    The solution to ODE #6 is: Eq(y(x), C1*exp(-x**2) + 1)


## Application

### Radioactive decay

**Slide 26**


```python
t = Symbol('t')  # Independent variable
N = Function('N')(t)  # Dependent variable N(t)
l = Symbol('l')  # Constant
# The ODE we wish to solve:
ode = diff(N,t) + l*N

print("ODE #7:") 
pprint(ode)
solution = dsolve(ode, N)
print("The solution to ODE #7 is:")
pprint(solution) 
```

    ODE #7:
             d       
    l⋅N(t) + ──(N(t))
             dt      
    The solution to ODE #7 is:
               -l⋅t
    N(t) = C₁⋅ℯ    


**Example:** 1 mole of carbon-14 at t=0


```python
l = 3.8394e-12
C = 6.02e23 * np.exp(l)

t_3d = np.arange(0, 1e12, 1e9)
n_3d = np.arange(0, 6.02e23, 6.02e20)
N, T = np.meshgrid(n_3d, t_3d)

dNdt = -l * N

t = np.linspace(0, 1e12, 1000)
n = C * np.exp(-l*t)
    
t_years = t/(3600*24*365.25)
   
fig = plt.figure(figsize=(16, 8))

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(T, N, dNdt, cmap='seismic', edgecolor='k', lw=0.25)
ax1.set_xlabel('t')
ax1.set_ylabel('N(t)')
ax1.set_zlabel('dN/dt')
ax1.set_title('Rate of change in number of C-14 nuclei', fontsize=14)

ax2 = fig.add_subplot(122)
ax2.plot(t_years, n, 'k', label='N(t) = Cexp(-lt)')
ax2.set_xlabel('Time (years)')
ax2.set_ylabel('Number of C-14 nuclei')
ax2.set_title("Number of C-14 nuclei over time", fontsize=14)
ax2.legend(loc='best', fontsize=12)
ax2.grid(True)

fig.tight_layout()
plt.show()
```

```{note}

The plane in the first graph shows that radioactive decay is independent of time, but only dependent on the number of radioactive nuclei present.

```

### Particle settling

**Slide 31**


```python
t = Symbol('t')  # Independent variable - time
v = Function('v')(t)  # Dependent variable v(t) - the particle velocity
# Physical constants
rho_f = Symbol('rho_f')  # Fluid density
rho_p = Symbol('rho_p')  # Particle density
eta = Symbol('eta')  # Viscosity
g = Symbol('g')  # Gravitational acceleration
a = Symbol('a')  # Particle radius

# The ODE we wish to solve.
ode = diff(v,t) - ((rho_p - rho_f)/rho_p)*g + (9*eta/(2*(a**2)*rho_p))*v

print("ODE #8:") 
pprint(ode)
solution = dsolve(ode, v)
print("The solution to ODE #8 is:")
pprint(solution) 
```

    ODE #8:
      g⋅(-ρ_f + ρₚ)   d          9⋅η⋅v(t)
    - ───────────── + ──(v(t)) + ────────
            ρₚ        dt            2    
                                 2⋅a ⋅ρₚ 
    The solution to ODE #8 is:
                                         ⎛       9⋅t  ⎞
                                       η⋅⎜C₁ - ───────⎟
                                         ⎜        2   ⎟
                2            2           ⎝     2⋅a ⋅ρₚ⎠
           - 2⋅a ⋅g⋅ρ_f + 2⋅a ⋅g⋅ρₚ + ℯ                
    v(t) = ────────────────────────────────────────────
                               9⋅η                     


**Example**: sand grain with density 2650kg/m3 and radius 1mm sinking in water.

Initial conditions: v=0 when t=0


```python
rho_f = 1000
rho_p = 2650
eta = 0.89
g = 9.81
a = 1e-3

C = -(2*a**2*rho_p)/(9*eta) * np.log((rho_p-rho_f)/rho_p)

v_ode = np.linspace(0, 10, 1000)
dvdt = (rho_p-rho_f)/rho_p - (9*eta*v_ode)/(2*a**2*rho_p)

t = np.arange(0, 0.005, 0.00005)
v = -2*a**2*g*rho_f + 2*a**2*g*rho_p + np.exp(eta*(C - 9*t/(2*a**2*rho_p)))/(9*eta)

fig = plt.figure(figsize=(16, 8))

ax1 = fig.add_subplot(121)
ax1.plot(v_ode, dvdt, 'r')
ax1.set_xlabel('v')
ax1.set_ylabel('dv/dt')
ax1.set_title("Plot of dv/dt vs v")
ax1.grid(True)

ax2 = fig.add_subplot(122)
ax2.plot(t, v, 'b')
ax2.set_xlabel('time (seconds)')
ax2.set_ylabel('velocity (m/s)')
ax2.set_title("Plot of velocity over time")
ax2.grid(True)

fig.tight_layout()
plt.show()
```


```python

```

537b16a205e911212ee8d9d94da8f052a730abb1
ipynb
Jupyter Notebook
notebooks/a_modules/math_methods_1/8_ODEs.ipynb
primer-computational-mathematics/book
305941b4f1fc4f15d472fd11f2c6e90741fb8b64
MIT
2020-08-02T07:32:14.000Z
2021-11-16T16:40:43.000Z
notebooks/a_modules/math_methods_1/8_ODEs.ipynb
primer-computational-mathematics/book
305941b4f1fc4f15d472fd11f2c6e90741fb8b64
MIT
2020-07-27T10:45:26.000Z
2020-08-12T15:09:14.000Z
notebooks/a_modules/math_methods_1/8_ODEs.ipynb
primer-computational-mathematics/book
305941b4f1fc4f15d472fd11f2c6e90741fb8b64
MIT
2020-08-05T13:57:32.000Z
2022-02-02T19:03:57.000Z
Qwen/Qwen-72B
1. YES
2. YES