import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import fsolve, curve_fit
from scipy.integrate import odeint, quad
from scipy.interpolate import interp1d
import sympy as sp
from IPython.display import display
sp.init_printing()
import pint; u = pint.UnitRegistry()
The following equation represents 1-D heat transfer in a heated rod of length $L=1$: $$\frac{d^2T}{dx^2} = -5432,$$ $$T(0) = 300,$$ $$T(L) = 300.$$
We can convert this second order ODE to two first order ODEs by making the substitution $H = dT/dx$. Then we have: $$\frac{dH}{dx} = -5432,$$
$$\frac{dT}{dx} = H,$$$$T(0) = 300,$$$$T(L) = 300.$$We can now solve these two coupled rate equations. But we only have an initial condition for $T$; we don't have an initial condition for H.
def rhsf(HT, x):
H = HT[0]
T = HT[1]
dHdx = -5432
dTdx = H
return np.array([dHdx, dTdx])
HT0 = np.array([3000, 300])
x = np.linspace(0, 1, 100)
HT = odeint(rhsf, HT0, x)
H = HT[:,0]
T = HT[:,1]
print(T[-1])
583.999999991825
def F(H0):
def rhsf(HT, x):
H = HT[0]
T = HT[1]
dHdx = -5432
dTdx = H
return np.array([dHdx, dTdx])
HT0 = np.array([H0[0], 300])
x = np.linspace(0, 1, 100)
HT = odeint(rhsf, HT0, x)
H = HT[:,0]
T = HT[:,1]
return T[-1] - 300
fsolve(F, 3000)
array([2716.00000001])