Скачать с ютуб How to: Solve an ODE in Python (Boundary Value Problem) в хорошем качестве

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