How to: Solve an ODE in Python (Boundary Value Problem)

Learn out to numerically solve an ordinary differential equation (ODE) in Python using a built in solver for boundary value problems: "scipy.integrate.solve_bvp()".

Here I discuss a 2nd order ODE with 2 boundary conditions (i.e. boundary value problem - BVP). A 1D steady-state heat transfer problem is used as an example, with various boundary conditions (fixed temperature, radiation) and with the option for heat generation. In order to solve the BPV, the 2nd order ODE must be converted into a system of 1st order ODEs in the format required by "scipy.integrate.solve_bvp()".

The Python script and PDF of the notes can be found here: https://www.hageslab.com/Resources.html#PythonAnchor

Here we are using "Spyder" IDE with the numpy, scipy, and matplotlib libraries

Script (for the radiation example):

import numpy as np #For basic math functions
import scipy.integrate as intg #For advanced math functions
import matplotlib.pyplot as plt #For plotting

#Define Constants
sigma = 5.67e-8 #[W/m^2/K^4]

#Set-up Paramters
qgen = 0 #[W/m**3]
k = 40 #[W/m/K]
T1 = 273.15 #[K]
epsilon = 0.8
Tsurr = 500+273.15 #[K]

#Set-up Grid
L = 1 #[m]
nodes = 100
x = np.linspace(0,L,nodes)

#Define Equations & Boundary Conditions
def f(x,y):
return np.vstack((y[1],np.full_like(x,-qgen/k)))

def bc(ya,yb):
return np.array([ya[0]-T1,k*yb[1]+sigma*epsilon*(yb[0]**4-Tsurr**4)])

#Inital Guess
y0 = np.zeros((2,x.size))

#Solve
sol = intg.solve_bvp(f,bc,x,y0)
T = sol.y[0]
dTdX = sol.y[1]

#Surface Temps
TS1 = T[0]
TS2 = T[-1]

#Compute Flux
q = -k*dTdX

#Plot
plt.figure(1,dpi=120)
plt.yscale('linear')
plt.xscale('linear')
#plt.xlim(0,1)
plt.ylim(273.15,500+273.15)
plt.title("Temperature Profile")
plt.xlabel("Distance / m")
plt.ylabel("T / K")
plt.plot(x,T)

plt.figure(2,dpi=120)
plt.yscale('linear')
plt.xscale('linear')
#plt.xlim(0,0.1)
#plt.ylim(0,100)
plt.title("Heat Flux")
plt.xlabel("Distance / m")
plt.ylabel("Flux / kW m$^{-1}$ K$^{-1}$")
plt.plot(x,q/1000)