The Scientific Method ===================== After some problems and examples, the method of separation of variables is introduced to apply the scientific method. Problems and Examples --------------------- We start with three different stories and two different examples. What they all have in common is that ordinary differential equations model those stories and examples. Problem 1: slowing down An attack submarine is cruising at 40 knots at a depth of 1000 feet when suddenly the reactor scrams. After 1 minute way has dropped to 30 knots. How long does the crew have to make repairs before forward motion falls below steerageway of 2 knots? Problem 2: cooling off A house furnace fails on a cold winter's evening when the outside (ambient) temperature of 20 degrees Fahrenheit. Although initially at 70 degree Fahrenheit, the inside temperature has fallen to 65 degrees Fahrenheit after 1 hour. How long before the inside temperature reaches the damaging temperature of 32 degrees Fahrenheit? Problem 3: population modeling Predict the future population of a developed country. Example 1: Planar Pendulum A mass :math:`m` is suspended on a rod of length :math:`L`. The motion of the pendulum stays in the plane, so we have only one parameter: the angle, measuring deviation from the stable equilibrium position. Example 2: Sliding Mass Spring The mass :math:`m` is sliding along a frictionless horizontal table restrained by a spring. The spring has constant :math:`k` and linear motion is assumed. We are interested in the displacement :math:`x(t)`, measured from the equilibrium, in function of time :math:`t`. Separation of Variables ----------------------- The solution to the Ordinary Differential Equation (ODE) .. math:: \frac{dy}{dx} = f(x) g(y), \quad y(x_0) = y_0 has a unique solution. We compute the solution via the separation of variables method. For example, consider: .. math:: \frac{dy}{dx} = \frac{x}{y}, \quad y(3) = 5. In applying ``SymPy.jl``, we first formulate the problem: :: x = Sym("x") y = SymFunction("y") ode = Eq(diff(y(x), x), x/y(x)) d x --(y(x)) = ---- dx y(x) To solve an ODE with ``SymPy``, we call ``dsolve``: :: sol = dsolve(ode) returns :: 2-element Vector{Sym}: Eq(y(x), -sqrt(C1 + x^2)), Eq(y(x), sqrt(C1 + x^2)) or equivalently :math:`y(x) = \pm \sqrt{C_1 + x^2}`. The parameter :math:`C_1` depends on the initial condition. Imposing the initial condition, the condition :math:`y(3) = 5` must be defined by a dictionary: :: initcond = Dict(y(3) => 5) Then we apply ``dsolve`` again on the ODE: :: solivp = dsolve(ode, ics=initcond) ________ / 2 y(x) = \/ x + 16 and now we see that :math:`y(x) = \sqrt{x^2 + 16}` is the unique solution. The scientific method consists in four steps: 1. Model the problem into an ODE. 2. Impose the data. 3. Solve the ODE. 4. Interpret the results. Often the steps may need to be repeated if the ODE is too simple to yield realistic results. Let us return to our first problem of slowing down: *An attack submarine is cruising at 40 knots at a depth of 1000 feet when suddenly the reactor scrams. After 1 minute way has dropped to 30 knots. How long does the crew have to make repairs before forward motion falls below steerageway of 2 knots?* :index:`Newton's law` equates the sum of the inertial forces to the sum of the external forces: * The inertial force is :math:`F = m a`, where :math:`m` is the mass and :math:`a` is the acceleration. * The external force is the water resistance :math:`R = -k v^2`, proportional to the square of the velocity, for some constant :math:`k`. So our model is .. math:: m \frac{d}{d~\!t} v(t) = - k ( v(t) )^2, where :math:`v(t)` is the velocity in function of time :math:`t`. The code to formulate the ODE with ``SymPy.jl`` is below: :: m, k, t = Sym("m, k, t") v = SymFunction(v) ode = Eq(m*diff(v(t),t), -k*v(t)^2) Let lump the constants :math:`m` and :math:`k` into one constant, dividing by :math:`m` and replacing :math:`k/m` by :math:`K`. :: lefteq = ode.lhs()/m K = Sym("K") righteq = subs(ode.rhs()/m, k/m=>K) newode = Eq(lefteq, righteq) d 2 --(v(t)) = -K*v (t) dt Before solving the ODE, we observe that our problem is an :index:`initial value problem`. Recall: *submarine is cruising at 40 knots* so the initial condition is :math:`v(0) = 40`. :: initcond = Dict(v(0)=>40) sol = dsolve(newode, ics=initcond) 1 v(t) = ---------- K*t + 1/40 What is :math:`K`? To answer this question, we impose the data. Recall: *after one minute way has dropped to 30 knots*, so we determine :math:`K` using :math:`v(1) = 30`. :: Kequ = subs(sol, t=>1, v(1) => 30) Ksol = solve(Kequ, K) The ``solve`` returns ``1-element Vector{Sym}`` so we use ``Ksol[1]`` in the next substitution: :: vsol = subs(sol, K=>Ksol[1]) 1 v(t) = -------- t 1 --- + -- 120 40 Interpreting the results, the original equation was *how long before the steerageway fall below 2 knots?* :: vequ = subs(vsol, v(t) => 2) solve(vequ, t) which returns 57. At 57 minutes the speed has dropped to 2 knots. So the crew has 56 minutes left. We continue to our second problem: *A house furnace fails on a cold winter's evening when the outside (ambient) temperature of 20 degrees Fahrenheit. Although initially at 70 degree Fahrenheit, the inside temperature has fallen to 65 degrees Fahrenheit after 1 hour. How long before the inside temperature reaches the damaging temperature of 32 degrees Fahrenheit?* :index:`Newton's cooling law` states that the rate at which the temperature changes is proportional to the gradient driving out the heat. Let :math:`T(t)` be the temperature at time :math:`t`. .. math:: \frac{d}{d~\!t} T(t) = k (T(t) - 20), \quad \mbox{for some constant } k. The left side of the ODE is the change in the temperature. The right side is the difference between the inside temperature and the outside temperature of 20 degree Fahrenheit. Formulating and solving the ODE with ``SymPy.jl`` goes as follows: :: T = SymFunction("T") ode = Eq(diff(T(t),t), k*(T(t) - 20)) The *initially at 70 degrees Fahrenheit* gives the initial condition: :: initcond = Dict(T(0) => 70) sol = dsolve(ode, ics=initcond) k*t T(t) = 50*e + 20 What is :math:`k`? Imposing the data, recall: *the inside temperature has fallen to 65 degrees Fahrenheit, after 1 hour.* :: kequ = subs(sol, t=>1, T(1)=>65) kval = solve(kequ, k) which gives :math:`\log(9/10)`, stored in an array. :: Tsol = subs(sol, k=>kval[1]) t*log(9/10) T(t) = 50*e + 20 Interpreting the results, we return to the original question: *How long before the inside temperature reaches the damaging temperature of 32 degrees Fahrenheit?* :: Tequ = subs(Tsol, T(t)=>32) tval = solve(Tequ, t) 1-element Vector{Sym}: log((6/25)^(1/log(9/10))) To interpret this exact result, we use a 64-bit floating-point approximation: :: Float64(tval[1]) 13.545077553292497 It takes about 13 and a half hour before it starts freezing inside. Our third problem asks to *Predict the future population of a developed country.* Considering the births, the growth of the population is proportional to its size. Let :math:`P(t)` denote the size of the population as function of time :math:`t`. .. math:: \frac{d}{d~\!t} P(t) = k P(t), \quad \mbox{for some constant } k. Formulating the ODE with ``SymPy.jl`` goes as follows: :: P = SymFunction("P") ode = Eq(diff(P(t),t), k*P(t)) dsolve(ode) k*t P(t) = C1*e The exponential is not realistic. Deaths occur. For a better model, consider the correction term, similar to the second term in a Taylor series. This second term is proportional to the square of the population size: :: logistics = Eq(diff(P(t),t), k*P(t) - K*(P(t))^2) This model is called the *logistics equation*. Solving the :index:`logistics equation`: :: Psol = dsolve(logistics) k*(C1 + t) k*e P(t) = ------------------- / k*(C1 + t) \ K*\e - 1/ To interpret this model, let us normalize the size of the population at the beginning of time to one. :: initcond = Dict(P(0)=>1) Psol1 = dsolve(logistics, ics=initcond) and then we obtain the expression: :: / / K \ \ | log|-----| | | \K - k/ | k*|t + ---------- | \ k / k*e P(t) = --------------------------- / / / K \ \ \ | | log|-----| | | | | \K - k/ | | | k*|t + ---------- | | | \ k / | K*\e - 1/ Applying ``simplify()``, we obtain .. math:: P(t) = \frac{k e^{k t}}{K e^{k t} - K + k}. In this expression, the two constants :math:`k` and :math:`K` are both positive. If we divide numerator and denominator by :math:`e^{k t}`: .. math:: P(t) = \frac{k}{K + (k-K) e^{-k t}}. As :math:`t \rightarrow \infty`, :math:`P(t) \rightarrow k/K`. This limit is called the *carrying capacity* of a society. Now we turn to the two mechanical models: a planar pendulum and a sliding mass connected to a spring, shown in :numref:`figPlanarPendulum` and in :numref:`figSlidingMass` respectively. .. _figPlanarPendulum: .. figure:: ./figPlanarPendulum.png :align: center A planar pendulum. We can model a planar pendulum in two ways: * Newton: tangential inertial force keeps the mass in place. The vertical component of the gravitational force is :math:`-m g \sin(\theta(t))`. .. math:: F = m L \frac{d^2}{d~\!t} \theta(t) = - m g \sin(\theta(t)). * Hamilton: kinetic and potential energy remain constant. The kinetic energy is :math:`\displaystyle \frac{1}{2} m L^2 \left( \frac{d}{d~\!t} \theta(t) \right)^2`. The gravitational force determines the potential energy is :math:`m g (L - L \cos(\theta(t)))`. Conservation of energy is expressed as .. math:: \frac{1}{2} m L^2 \left( \frac{d}{d~\!t} \theta(t) \right)^2 + m g (L - L \cos(\theta(t))) = E, for some constant :math:`E`. Consider a mass :math:`m` sliding along a frictionless horizontal table, restrained by a spring. The displacement from equilibrium is :math:`x(t)`, as shown in :numref:`figSlidingMass`. .. _figSlidingMass: .. figure:: ./figSlidingMass.png :align: center A sliding mass connected to a spring. By :index:`Newton's law`, the inertial force is :math:`\displaystyle m \frac{d^2}{d~\!t} x(t)`. The restoring force of the spring is proportional to :math:`x(t)`, so we obtain .. math:: m \frac{d^2}{d~\!t} x(t) = -k x(t), \quad \mbox{for some constant } k. Connecting several masses with springs serves as a good model for car traffic, a rich source of industrial problems. Proposal for a Project Topic ---------------------------- How well does the logistic model fit the population of the U.S.A.? * Gather census data of the U.S. population since 1945. * Fit the constants :math:`k` and :math:`K`. * What is the carrying capacity? Interpret the results. Exercises --------- 1. A coasting rowboat experiences * subsurface resistance proportional to :math:`v^2`, and * bow-wave surface resistance proportional to :math:`v^3`. Model the velocity :math:`v` of this craft. Will it ever come to rest? 2. A can of pop at :math:`34^{\circ}\mbox{F}` is removed from a refrigerator, placed on a counter for awhile in a room at temperature of :math:`75^{\circ}\mbox{F}`. After 2 minutes the temperature of the can has risen to :math:`35^{\circ}\mbox{F}`. Predict the future evolution of the temperature of the can. When will the temperature of the can reach :math:`70^{\circ}\mbox{F}`? 3. Show that the equation .. math:: \frac{1}{2} m L^2 \left( \frac{d}{d~\!t} \theta(t) \right)^2 + m g (L - L \cos(\theta(t))) = E is equivalent to the second order differential equation derived with Newton's rule. *Hint:* apply :math:`\displaystyle \frac{d}{d~\!t}` to the equation and factor.