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

Veröffentlicht am: 07 September 2022
auf dem Kanal: 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)

Auf dieser Seite können Sie das Online-Video How to: Solve an ODE in Python (Boundary Value Problem) mit der Dauer stunde minuten sekunde in guter Qualität ansehen, das der Benutzer HagesLab 07 September 2022 hochgeladen hat, den Link mit Freunden und Bekannten teilen, dieses Video wurde auf Youtube bereits 10,363 Mal angesehen und es wurde von 211 den Zuschauern gefallen. Viel Spaß beim Betrachtenden Zuschauern gefallen!

291