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

Publicado el: 07 septiembre 2022
en el canal de: HagesLab

10,363

211

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.ht...

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)

En esta página del sitio puede ver el video en línea How to: Solve an ODE in Python (Boundary Value Problem) de Duración hora minuto segunda en buena calidad , que subió el usuario HagesLab 07 septiembre 2022, comparta el enlace con amigos y conocidos, en youtube este video ya ha sido visto 10,363 veces y le gustó 211 a los espectadores. Disfruta viendo!

291