In [2]:
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()

Sample problem

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.

Part a

  • Find T(x) and H(x) by solving the two rate equations assuming $H(0) \equiv H0= 3000$.
  • (Note, because we had to guess H(0), $T(L)\ne 300$ like it should be. We'll fix this in Part b.)
In [4]:
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

Part b

  • We can treat all of part A as a solution of one "nonlinear equation" in one unknown.
    • The unknown is H0.
    • The "nonlinear equation" is not some formula we write down as usual.
      • Rather, it is a function: put in a value of H0, solve the ODEs for T(x) and H(x), evaluate and return T(L).
      • Call this function f(H0). You put in H0, and you compute and returnn T(L).
    • So, we are solving for f(H0) = 300.
    • Or, F(H0) = f(H0)-300 = 0.
  • Wrap all of Part a in a function F(H0), and solve F(H0)=0 for H0.
  • Report the value rounded as ####
In [5]:
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)
    
    
Out[5]:
array([2716.00000001])
In [ ]: