In [2]:
import numpy             as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize      import fsolve, curve_fit
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 [ ]: