.. _ode: Numerical Differential Equations ================================ .. include:: ../replace.txt .. index:: differential equations, ode .. topic:: Reading Assignment Please read chapter 8 of *Physical Modeling in MATLAB*, by Allen B. Downey [DOWNEY11]_. .. index:: Numerical Differential Equations, ODE An ordinary differential equation (ODE) is an equation of the form .. math:: y'(t) = \frac{dy}{dt} = f(t, y), \text{ with } y(t_0) = y_0. We are given the rate of change (the derivative function :math:`f(t, y)`) of :math:`y(t)` from which we solve for :math:`y(t)`. Because we are also given the initial value :math:`y(t_0) = y_0`, this class of ODE problem is called an initial value problem (IVP). We will consider cases of the function :math:`y(t)` that represents a single equation and also a system of equations. Differential equations are important to engineers and scientists because changes to the state of dynamic systems are described by differential equations. Examples of these dynamic systems include objects that move; wires carrying fluctuating electrical current; reservoirs of fluid with measurable pressure, temperature, and volume; or anything where the state of the system can be quantified numerically. Numerical methods are especially needed to solve higher order and nonlinear ODEs. As is the case with integration, our ability to properly evaluate some ODEs using analytic methods is less than satisfying [STRANG14]_. So numerical methods for solving ODEs is an important tool in the engineer's computational toolbox. This is our third encounter with ODEs. In :ref:`ballup`, we used |M|'s Symbolic Math Toolbox to find the algebraic equations from the differential equations associated with the effect of gravity on a ball thrown upwards. Then in :ref:`eigODE` we considered systems of first order linear ODEs that may be solved analytically with eigenvalues and eigenvectors. Now, we will consider numerical methods to find and plot the solutions to ODEs. We will cover the basic algorithmic concepts of how numerical ODE solvers work and review the efficient suite of ODE solver provided by |M|. Recall the definition of a derivative: .. math:: \frac{dy}{dt} = \lim_{h \rightarrow 0} \frac{y(t + h) - y(t)}{h}. Considered in terms of discrete samples of :math:`y`, the derivative becomes .. math:: y_k' = \frac{y_{k+1} - y_k}{h}, where :math:`h` is the distance between sample points, :math:`h = t_{k+1} - t_k`. Rearranged to find :math:`y_{k+1}` from the previous value, :math:`y_k`, and the slope of the function at :math:`t_k`, we have .. math:: y_{k+1} = y_k + h\,y'_k :label: eq-euler1 This equation is the basis for the simplest ODE solver, Euler's method. However, due to concerns about excessive computation and round-off errors, we are limited on how small :math:`h` can be. So the accuracy of :eq:`eq-euler1` is not acceptable for most applications. Since differential equations are given in terms of an equation for the derivative of :math:`y`, :math:`y' = f(t, y)`, and the initial value, :math:`y_0`, the general discrete solution comes from the fundamental theorem of calculus, which applies the antiderivative of a function to a definite integral. .. math:: y_{k+1} = y_k + \int_{t_k}^{t_{k+1}} f(t, y) dt :label: eq-genODE Then the numerical approximation comes from the algorithm used to calculate the definite integral. .. math:: I_k = \int_{t_k}^{t_{k+1}} f(t, y) dt :label: eq-odeInt At each step, the ODE solver uses a straight line with a slope of :math:`I_k/h` to advance :math:`y_k` to :math:`y_{k+1}`. .. _euler: Euler's Method --------------- Euler's method uses a Reimann integral to approximate :eq:`eq-odeInt` with constant height rectangular subintervals. It is a reasonable solution if :math:`f(t, y)` is constant or nearly constant. However, the results from Euler's methods applied to most ODEs are not accurate enough for most applications unless the step size is small enough relative to the rate of change of the :math:`f(t, y)` function. .. \footnote{This is the \emph{forward} Euler method. A variant to this algorithm is the implicit backward Euler method, which we will use in \cref{backEuler}.} .. math:: \int_{t_k}^{t_{k+1}} f(t, y) dt \approx h\,f(t_k, y_k) .. math:: y_{k+1} = y_k + h\,f(t_k, y_k) Euler's method is a perfect match if we let :math:`f(t, y) = a`, then .. math:: y_{k+1} = y_k + h\,a. Then by induction, .. math:: \spalignsys{y_1 = {y_0 + h\,a} {}; y_2 = {y_1 + h\,a} {= y_0 + 2\,h\,a}; \. \vdots; y_k = {y_0 + k\,h\,a} {}}. The ``eulerODE`` function listed below is an implementation of Euler's method. The first three arguments passed to ``eulerODE`` are the same as is used with the |M| ODE solvers. An additional argument for the number of values to find is also passed to ``eulerODE``. First, define a function or function handle. The function defines the derivative of the solution. It should take two arguments, :math:`t` and :math:`y`. Although in many cases, the function only operates on one variable---usually the second, ``y``, variable. The function should accept and return vector variables. .. math:: \frac{dy}{dt} = f(t, y) Next, specify a ``tspan`` argument, which is a row vector of the initial and final :math:`t` values (``[t_0 t_final]``). Then, the value of :math:`y` at the initial time is specified :math:`y(t_0)`. The final argument is :math:`n`, the number of values to find. The returned data from all of the ODE solvers are the :math:`t` evaluation points and the :math:`y` solution. .. \footnote{We use $t$ for time as a variable of the ODE function, but a spatial variable, such as $x$, could also be used.} :: >> f = @(t, y) 5 >> [t, y] = eulerODE(f, [0 10], 5, 4) t = 0 3.3333 6.6667 10.0000 y = 5.0000 21.6667 38.3333 55.0000 :numref:`fig-ODEcompare` shows a plot comparing the relative accuracy of Euler's method to other fixed step size algorithms. :: function [t, y] = eulerODE(f, tspan, y0, n) % EULERODE - ODE solver using Euler's method. % % f - vectorized derivative function, y' = f(t, y) % tspan - row vector [t0 tf] for initial and final time % y0 - initial value of y % n - number of y values to find h = (tspan(2) - tspan(1))/(n-1); t = linspace(tspan(1), tspan(2), n); y = zeros(1, n); y(1) = y0; for k = 1:n-1 y(k+1) = y(k) + h*f(t(k), y(k)); end .. \label{code-eulerODE} .. \cref{code-eulerODE} Heun's Method --------------- Heun's method uses trapezoid rule integration to approximate :eq:`eq-odeInt`. .. math:: :label: eq-Heun1 y_{k+1} = y_k + \frac{h}{2}\left[f(t_k, y_k) + f(t_{k+1}, y_{k+1})\right] However, :eq:`eq-Heun1` makes use of :math:`y_{k+1}`, which is not yet know. So Euler's method is first used to make an initial estimate of :math:`y_{k+1}`. Because the algorithm makes an initial prediction of :math:`y_{k+1}` and then updates the estimate with a later equation, Heun's method belongs to a class of algorithms called *predictor--corrector methods*. .. math:: y_{k+1} = y_k + h\,f(t_k, y_k) The complete Heun method approximation equation is given by :eq:`eq-Heun2`. .. math:: :label: eq-Heun2 y_{k+1} = y_k + \frac{h}{2}\left[f(t_k, y_k) + f(t_{k+1}, (y_k + h\,f(t_k, y_k))\right] The ``heunODE`` function listed below is called with the same arguments as the ``eulerODE`` function. :: function [t, y] = heunODE(f, tspan, y0, n) % HEUNODE - ODE solver using Heun's method. % % f - vectorized derivative function, y' = f(t, y) % tspan - row vector [t0 tf] for initial and final time % y0 - initial value of y % n - number of y values to find h = (tspan(2) - tspan(1))/(n-1); t = linspace(tspan(1), tspan(2), n); y = zeros(1, n); y(1) = y0; for k = 1:n-1 fk = f(t(k), y(k)); yk1 = y(k) + h*fk; y(k+1) = y(k) + (h/2)*(fk + f(t(k+1), yk1)); end .. \label{code-heunODE} Runge-Kutta Method -------------------- The Runge-Kutta method uses Simpson's rule integration to approximate :eq:`eq-odeInt`. There are actually several algorithms that are considered to be part of the Runge-Kutta family of algorithms. The so called *classic* algorithm used here is an order--4 algorithm with a comparatively simple expression. Runge-Kutta algorithms are weighted sums of coefficients calculated from evaluations of :math:`f(t, y)` (the derivative function). Later coefficients make use of previous coefficients, so the coefficients at each evaluation point are first calculated sequentially before the solution value can be determined. Runge-Kutta algorithms are the most accurate of the fixed step size algorithms. They are only surpassed by adaptive step size algorithms, which typically make use of Runge-Kutta equations. .. math:: \begin{array}{rl} s_1 &= f(t_k, y_k) \\ \hfill s_2 &= f(t_k + \frac{h}{2}, y_k + \frac{h}{2}\,s_1) \\ \hfill s_3 &= f(t_k + \frac{h}{2}, y_k + \frac{h}{2}\,s_2) \\ \hfill s_4 &= f(t_k + h, y_k + h\,s_3) \\ \hfill y_{k + 1} &= y_k + \frac{h}{6}\, \left( s_1 + 2s_2 + 2s_3 + s_4 \right) \end{array} The ``rk4ODE`` function listed below is called with the same arguments as the ``eulerODE`` and ``heunODE`` functions. :: function [t, y] = rk4ODE(f, tspan, y0, n) % RK4ODE - ODE solver using the classic Runge-Kutta method. % % f - vectorized derivative function, y' = f(t, y) % tspan - row vector [t0 tf] for initial and final time % y0 - initial value of y % n - number of y values to find h = (tspan(2) - tspan(1))/(n - 1); t = linspace(tspan(1), tspan(2), n); y = zeros(1, n); y(1) = y0; for k = 1:n-1 s1 = f(t(k), y(k)); s2 = f(t(k)+h/2, y(k)+s1*h/2); s3 = f(t(k)+h/2, y(k)+s2*h/2); s4 = f(t(k+1), y(k)+s3*h/2); y(k+1) = y(k) + (h/6)*(s1 + 2*s2 + 2*s3 + s4); end .. _ode-alg-comp: Algorithm Comparison ----------------------- We will use the well known decaying exponential ODE to compare the three methods that use a fixed step size. We will also reference this simple ODE when considering stability constraints on the step size in :ref:`ode-stability`. The ODE .. math:: y' = \frac{dy(t)}{dt} = -a\,y,\; \text{where}\; a > 0,\: \text{and}\: y(0) = c has the solution :math:`y(t) = c\,e^{-at}`. For our comparison, :math:`a = 1`, and :math:`y(0) = 1`. We can see in :numref:`fig-ODEcompare` that, as expected, the classic Runge-Kutta method preforms the best followed by Heun's method and then Euler's method. :: >> f = @(t, y) -y; >> [t, y_euler] = eulerODE(f, [0 4], 1, 8); >> y_exact = exp(-t); >> [t, y_heun] = heunODE(f, [0 4], 1, 8); >> [t, y_rk4] = rk4ODE(f, [0 4], 1, 8); .. _fig-ODEcompare: .. figure:: ODEcompare.png :width: 80% :align: center Comparison of fixed interval ODE solvers. .. 9 cm wide .. note:: All of the algorithms shown in :numref:`fig-ODEcompare` have very good performance if :math:`h` is significantly reduced. Finding noticeable differences in the plot is only achieved by using a larger step size. .. _ode-stability: Stability Constraints ------------------------ Stability and stability constraints are important to us because we do not want our numerical solution to give unstable and wrong answers. A user without knowledge of the stability constraints could specify a step size for the fixed step size algorithm that is too large to get a stable numerical solution. The stability constraint is especially concerning with stiff systems where one term in the solution may place a much higher constraint on the step size than other terms in the solution. .. \footnote{Stiff ODEs are addressed in \cref{stiffODE}.} We return to the Euler method and an exponentially decaying ODE to find a stability constraint related to :math:`h`, the distance between consecutive evaluation points of ODE solvers [SHEN15]_. With :math:`y' = -a\,y`, where :math:`a > 0`, :math:`y(0) = 1`, the exact solution is :math:`y(t) = e^{-a\,t}`. A property of :math:`y(t)` that needs to be preserved is that :math:`y \to 0` as :math:`t \to \infty`. Referring to :ref:`euler`, The forward Euler solution is :math:`y_{k+1} = y_k - h\,a\,y_k`, or :math:`y_{k+1} = (1 - ha)\,y_k`. By induction, .. math:: y_k = (1 - ha)^k\,y_0 A stability constraint is needed to ensure that the property of :math:`y \to 0` as :math:`k \to \infty` is preserved. .. math:: \begin{array}{rl} &\abs{1 - h\,a} < 1 \\ -1 &< (1 - h\,a) < 1 \\ -2 &< -h\,a < 0 \\ &h < \frac{2}{a} \end{array} In the following example :math:`2/a = 2/0.5 = 4`, but :math:`h = 5` and the results are unstable. :: >> fexp = @(t, y) -0.5*y; >> [t, y_unstable] = eulerODE(fexp, [0 20], 1, 5) t = 0 5 10 15 20 y_unstable = 1.0000 -1.5000 2.2500 -3.3750 5.0625 The other ODE solvers evaluate the function at points between the set evaluation points, so the stability constraint for Euler's method is a sufficient constraint for them. .. _backEuler: Implicit Backward Euler Method -------------------------------- Implicit ODE solvers look at :eq:`eq-genODE` from the perspective of finding :math:`y_k` from the yet unknown :math:`y_{k+1}`. So the terms of :eq:`eq-genODE` are reordered and the limits of the integral are reversed. The resulting algorithm removes the stability constraint on the step size. .. math:: y_k = y_{k+1} + \int_{t_{k+1}}^{t_k} f(t, y) dt :label: eq-implicitODE Again, we will use the ODE :math:`y' = -a\,y`, where :math:`a > 0`, and :math:`y(0) = 1`. The exact solution is :math:`y(t) = e^{-a\,t}`. As before, we want to find any stability constraints that are required to ensure that :math:`y \to 0` as :math:`k \to \infty`. We can solve the decaying exponential ODE problem numerically by applying Euler's method to :eq:`eq-implicitODE`. .. math:: :label: eq-implicitEuler \begin{array}{rl} y_k &= y_{k+1} + h\,a\,y_{k+1} \\ y_k &= (1 + h\,a)\,y_{k+1} \\ y_{k+1} &= \frac{1}{1 + h\,a}\,y_k \\ y_{k+1} &= \left(\frac{1}{1 + h\,a}\right)^k y_0 \\ \end{array} To satisfy the stability constraint, we require that .. math:: \left\vert\frac{1}{1 + h\,a}\right\vert < 1, which always holds because both :math:`h` and :math:`a` are positive numbers. So the implicit backward Euler method is unconditionally stable. .. _adaptiveODE: Adaptive Algorithms -------------------- .. index:: adaptive ODE algorithms Numerical ODE solvers need to use appropriate step sizes between evaluation points to achieve high levels of accuracy while also quickly finding solutions across the evaluation interval. Adaptive ODE solvers reduce the step size, :math:`h`, where needed to meet the accuracy requirements. To improve performance, they also increase the step size when the accuracy requirements can be satisfied with a larger step size. One of two strategies can be used to determine the step size adjustments. As was done with adaptive numerical methods for computing integrals in :ref:`adaptIntegral`, the solution from two half steps could be compared to the solution from a full step. Based on the comparison, the step size can be reduced, left the same, or increased. However, points of an ODE solution are found sequentially, so the strategy of splitting a wide region into successively smaller regions would be computationally slow. Instead, the current step size is used to find the next solution (:math:`y_{k+1}`) using two algorithms where one algorithm is regarded as more accurate than the other. The guiding principle is that if the solutions computed by the two algorithms are close to the same, then the step size is small enough or could even be increased. But if the difference between the solutions is larger than a threshold, then the step size should be reduced. If the solutions differ significantly, then the current evaluations are recalculated with a smaller step size. The algorithm advances to the next evaluation point only when the differences between the solutions is less than a specified tolerance. .. _rkf45: The RKF45 Method ^^^^^^^^^^^^^^^^^^^^ .. index:: RKF45 method The Runge-Kutta-Fehlberg method (denoted as RKF45) uses fourth and fifth order Runge-Kutta ODE solvers to adjust the step size at each evaluation point. The fourth order equation is more complex than is used by the fourth order classic Runge-Kutta algorithm. But the computed coefficients for the fourth order equation are also used for the fifth order equation, which only requires one additional coefficient. The :math:`s_2` coefficient is used in finding later coefficients, but is not part of either the fourth or fifth order solution equations. In the ``rkf45ODE`` function, the fourth order solution values (RK4) are stored in the ``z`` array, and the fifth order solution values (RKF45) are stored in the ``y`` array, which is returned from the function. The coefficients and solution values are calculated using the following equations. .. math:: \begin{array}{rl} &s_1 = h\,f\left(t_k, y_k\right) \\ \hfill &s_2 = h\,f\left(t_k + \frac{1}{4}h, y_k + \frac{1}{4}s_1\right) \\ \hfill &s_3 = h\,f\left(t_k + \frac{3}{8}h, y_k + \frac{3}{32}s_1 + \frac{9}{32}s_2\right) \\ \hfill &s_4 = h\,f\left(t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}s_1 - \frac{7200}{2197}s_2 + \frac{7296}{2197}s_3\right) \\ \hfill &s_5 = h\,f\left(t_k + h, y_k + \frac{439}{216}s_1 - 8s_2 + \frac{3680}{513}s_3 - \frac{845}{4104}s_4\right) \\ \hfill &s_6 = h\,f\left(t_k + \frac{1}{2}h, y_k - \frac{8}{27}s_1 + 2s_2 - \frac{3544}{2565}s_3 + \frac{1859}{4104}s_4 - \frac{11}{40}s_5\right) \\ \hfill &z_{k+1} = y_k + \frac{25}{216}s_1 + + \frac{1408}{2565}s_3 + \frac{2197}{4101}s_4 - \frac{1}{5}s_5 \\ \hfill &y_{k+1} = y_k + \frac{16}{135}s_1 + \frac{6656}{12825}s_3 + \frac{28561}{56430}s_4 - \frac{9}{50}s_5 + \frac{2}{55}s_6 \end{array} There are variances in the literature on the implementation details of the RKF45 method. The ``rkf45ODE`` function listed below is adapted from the description of the algorithm in the numerical methods text by Mathews and Fink [MATHEWS05]_. Some adaptive ODE algorithms make course changes to the step size such as either cutting it in half or doubling it. The ``rkf45ODE`` function first finds a scalar value, :math:`s`, from a tolerance value and the difference between the solutions from the two algorithms. The new step size comes from multiplying the current step size by the scalar. The scalar will be one if the difference between the solutions is half of the tolerance. The algorithm advances to the next step if the difference between the two solutions is less than the tolerance value, which occurs when the scalar is greater than 0.84. Otherwise, the current evaluation is recalculated with a smaller step size. The multiplier is restricted to not be more than two. So the step size will never increase to more than twice its current size. The tolerance and the threshold for advancing to the next point are tunable parameters. .. math:: s = \left(\frac{tolerance}{2\,\norm{z_{k+1}- y_{k+1}}} \right)^{\frac{1}{4}} The initial value (:math:`y_0`) passed to the ``rkf45ODE`` function can be either a scalar for a single equation ODE or a column vector for a system of ODEs. Since the number of steps to be taken is not known in advance, the initial step size is a function argument instead of the number of solution values to find. :: function [t, y] = rkf45ODE(f, tspan, y0, h0) % rk5ODE - Adaptive interval ODE solver using the % Runge-Kutta-Fehlberg method. % % f - vectorized derivative function, y' = f(t, y) % tspan - row vector [t0 tf] for initial and final time % y0 - initial value of y, scalar or column vector % h0 - intial step size tol = 2e-5; h = h0; nEst = ceil((tspan(2) - tspan(1))/h); m = size(y0, 1); t = zeros(1, 2*nEst); % estimated preallocation y = zeros(m, 2*nEst); z = zeros(m, 2*nEst); t(1) = tspan(1); y(:,1) = y0; % RKF45 results z(:,1) = y0; % RK4 results k = 1; while t(k) < tspan(2) s1 = h*f(t(k), y(:,k)); s2 = h*f(t(k) + h/4, y(:,k) + s1/4); s3 = h*f(t(k) + h*3/8, y(:,k) + s1*3/32 + s2*9/32); s4 = h*f(t(k) + h*12/13, y(:,k) + s1*1932/2197 ... - s2*7200/2197 + s3*7296/2197); s5 = h*f(t(k) + h, y(:,k) + s1*439/216 - s2*8 ... + s3*3680/513 - s4*845/4104); s6 = h*f(t(k) + h/2, y(:,k) - s1*8/27 + s2*2 ... - s3*3544/2565 + s4*1859/4104 - s5*11/40); z(:,k+1) = y(:,k) + s1*25/216 + s3*1408/2565 ... + s4*2197/4101 - s5/5; % RK4 y(:,k+1) = y(:,k) + s1*16/135 + s3*6656/12825 ... + s4*28561/56430 - s5*9/50 + s6*2/55; % RKF5 % no divide by 0, h at most doubles algDiff = max(tol/32, norm(z(:,k+1) - y(:,k+1))); s = (tol/(2*algDiff))^0.25; h = h*s; if s > 0.84 t(k+1) = t(k) + h; k = k + 1; end end % delete unused preallocated memory t(k+1:end) = []; y(:,k+1:end) = []; :numref:`fig-compRk4Rkf45` shows that the adaptive Runge-Kutta-Fehlberg (RKF45) method is more accurate than the classic Runge-Kutta method (RK4) with the initial step size of RKF45 being the same as the fixed step size used by RK4. For most ODEs, the ability of RKF45 to reduce the step size where needed likely contributes more to its improved accuracy than can be attributed to using a higher order equation. :: >> f = @(t,y)-y; >> [t, y_rkf45] = rkf45ODE(f, [0 4], 1, 0.5); >> y_exact = exp(-t); >> hold on >> plot(t, y_exact) >> plot(t, y_rkf45, ':') >> [t, y_rk4] = rk4ODE(f, [0 4], 1, 8); >> plot(t, y_rk4, '--') >> hold off, axis tight, xlabel('t') >> ylabel('y(t)'), legend('exact', 'RKF45', 'RK4') .. _fig-compRk4Rkf45: .. figure:: compRk4Rkf45.png :align: center :width: 80% The adjustable step size of the RKF45 algorithm gives accurate results with good performance. .. 9 cm wide MATLAB's ``ode45`` ^^^^^^^^^^^^^^^^^^^ Although the adaptive ``rkf45ODE`` function has good results, it is not as accurate as the |M| ``ode45`` function. Using the same comparison test shown in :numref:`fig-compRk4Rkf45`, the ``ode45`` function had a maximum error of :math:`6{\times}10^{-6}`, and a plot of its solution appears to be directly on top of the exact solution plot. |M|'s ``ode45`` is the ODE solver that should be tried first for most problems. It has good accuracy and performance for non-stiff ODEs. Like the RKF45 method, it uses an adaptive step size algorithm based on the fourth and fifth order Runge-Kutta formula. .. \footnote{Stiff ODEs are addressed in \cref{stiffODE}.} The ``ode45`` function is called with the same arguments as the ODE solvers previously described, except neither the initial step size nor the number of evaluation points is given. The |M| ODE solvers also offer additional flexibility regarding the evaluation points. If the ``tspan`` argument contains more than two values, then ``tspan`` is taken to be a vector listing the evaluation points to use. As before, we often save two return values from ``ode45``, ``[t, y] = ode45(f, tspan, y0)``. The return values in are column vectors rather than row vectors as was the case with the simple functions previously used. If no return values are specified, ``ode45`` will make a plot of the :math:`y` values against the :math:`t`--axis. The |M| ODE solvers have several additional options and tunable parameters. Rather than passing those individually to the solver, a data structure holding the options is passed as an optional argument. The data structure is created, and modified using the ``odeset`` command with name--value pairs. Refer to the |M| documentation for ``odeset`` for a full list of the options and their descriptions. :: >> options = odeset('RelTol', 1e-5); % create >> options = odeset(options, 'AbsTol',1e-7); % modify >> [t, y] = ode45(f, tspan, y0, options); % use .. _logisticEQ: The Logistic Equation ---------------------- This example problem is a model of growth slowed down by competition. The competition is within the population for limited resources. The growth gives a :math:`y` term in the ODE---the growth is proportional to the population size. The competition (:math:`y` against :math:`y`) adds a :math:`-y^2` term. .. math:: \frac{dy}{dt} = a\,y - b\,y^2. The solution is modeled with positive and negative values for :math:`t`. The analytic solution is: ([STRANG14]_ page 55) .. math:: y(t) = \frac{a}{b + d\,e^{-t}}. Where, :math:`d = \frac{a}{y(0)} - b`. In the distant past, (:math:`t = -\infty`), the population size is modeled as zero. The steady state (:math:`t = \infty`) population size is :math:`y(t) = K = \frac{a}{b}`, which is the sustainable population size. The half way point, can be modeled as :math:`y(0) = a/2b`. We will compute the positive half of the solution and then reflect the values to get the negative time values. :: % File: logistic_diffEq.m % Logistic Equation: a = 1, b = 1 f = @(t, y) y - y.^2; tspan = [0 5]; y0 = 1/2; [t, y] = ode45(f, tspan, y0); figure, plot(t,y) t1 = flip(-t); y1 = flip(1-y); t2 = [t1;t]; y2 = [y1;y]; analytic = 1./(1 + exp(-t2)); figure, hold on plot(t2, y2, 'o') plot(t2, analytic) hold off, title('Logistic Equation') xlabel('t'), ylabel('y') legend('ode45 Numeric Solution','Analytic Solution', ... 'Location', 'northwest'); .. _fig-logistic: .. figure:: logistic.png :align: center :width: 80% Numerical and analytical solution to the logistic ODE. .. _stiffODE: Stiff ODEs ---------------- A problem is stiff if the solution contains multiple terms with different step size requirements. A slower moving component may dominate the solution; however, a less significant component of the solution may change much faster. Adaptive step size algorithms must take small steps to maintain accurate results. Because stiff ODEs problems have solutions with multiple terms, they come from either second or higher order ODEs or from systems of ODEs. We will now look at an example of a moderately stiff system of ODE equations. Then we will see how implicit solvers solve stiff systems of ODE and are unconditionally stable. Then we will loot at |M|'s stiff ODE solvers, such as ``ode15s``. .. _stiffSys: Stiff ODE Systems ^^^^^^^^^^^^^^^^^^ Here we show a system of ODEs that have a moderately stiff solution. The reader is referred to :ref:`eigODE` for a review of systems of first order, linear ODEs. .. math:: \spalignsys{x' = -25\,x + 24\,y{,} {\qquad} x(0) = 2; y' = 24\,x - 25\,y{,} {\qquad} y(0) = 1} It is more convenient to express systems of ODEs using vectors and matrices. .. math:: \bm{u}(t) = \vector{x(t); y(t)} \qquad \mathbf{A} = \mat{-25 24; 24 -25} \qquad \bm{u'} = \mathbf{A}\,\bm{u} \qquad \bm{u}_0 = \vector{1; 2} The eigenvectors, eigenvalues, and :math:`c` coefficients are found with some help from |M|. .. math:: \bm{x}_1 = \frac{1}{\sqrt{2}}\vector{1; -1} \quad \bm{x}_2 = \frac{1}{\sqrt{2}}\vector{1; 1} \quad \spaligntabular{l}{{$\lambda_1 = -49$}; {$\lambda_2 = -1$}} \quad \spaligntabular{l}{{$c_1 = -\sqrt{2}\cdot 0.5$}; {$c_2 = \sqrt{2}\cdot 1.5$}} .. refer eq from eigODE The solution to the system of ODEs then comes from the general solution for systems of first order, linear ODEs as given in :ref:`eigODE`. .. math:: \bm{u}(t) = \spalignsys{x(t) = -0.5\,e^{-49\,t} + 1.5\,e^{-t}; y(t) = 0.5\,e^{-49\,t} + 1.5\,e^{-t}} This system is stiff because the exact solution has terms with largely different exponent values coming from the eigenvalues of :math:`\bf{A}`. The stability constraints on the :math:`e^{-49\,t}` terms require :math:`h < 2/49 \approx 0.0408`. From the beginning of the evaluation, those terms will move towards zero much faster than the :math:`e^{-t}` terms. We can see from :numref:`fig-stiffSys` and :numref:`fig-stiffError` that an adaptive algorithm such as RKF45 or |M|'s ``ode45`` will significantly reduce the step size from the beginning of the evaluations until the :math:`e^{-49\,t}` term has sufficiently gone to zero. Using a longer evaluation time, the errors between the RKF45 solution and the exact solution go to near zero and the step sizes get much larger. :: >> A = [-25 24; 24 -25]; >> fA = @(t, u) A*u; >> u0 = [1;2]; >> [t, u] = rkf45ODE(fA, [0 1], u0, 0.04); >> x = @(t) -0.5*exp(-49*t) + 1.5*exp(-t); >> y = @(t) 0.5*exp(-49*t) + 1.5*exp(-t); >> figure, hold on >> plot(t, u(1,:), 'o') >> plot(t, u(2,:), '+') >> plot(t, x(t), ':') >> plot(t, y(t), '--') >> hold off >> axis tight >> legend('u(1,:)', 'u(2,:)', 'x(t)', 'y(t)') >> xlabel('t') .. sideByside .. _fig-stiffSys: .. figure:: stiffSys.png :align: center :width: 80% Stiff ODE system---the :math:`\pm 0.5\,e^{-49\,t}` terms quickly goes to zero and the :math:`x(t)` and :math:`y(t)` solutions approximately merge to :math:`1.5\,e^{-t}`. :: >> U = [x(t); y(t)]; >> figure; >> plot(t, vecnorm(U - u)) >> hold on >> plot(t(1:end-1), diff(t), 'o:') >> axis tight >> xlabel('t') >> legend('error vector', 'h', 'Location', 'southeast') >> hold off .. _fig-stiffError: .. figure:: stiffError.png :align: center :width: 80% Error and step size for stiff ODE system. First curve---the norm of the error vector between the exact solution and that found from the RKF45 solver. Second curve---the step size, :math:`h`, for the stiff ODE system. .. _IBEsys: Implicit Solvers for Systems of ODEs ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The implicit backward Euler method is used here to solve stiff systems of ODEs. The initial discussion parallels that of :ref:`backEuler` but the equations are extended to systems of ODEs. As with the implicit solver for a single ODE equation, the implicit backward Euler (IBE) method for systems of ODEs is also unconditionally stable. We will again use a decaying exponential ODE such as described in :ref:`stiffSys`. As was the case in :ref:`stiffSys`, the :math:`\bf{A}` matrix is required be a symmetric matrix with real, negative eigenvalues. .. math:: \bm{u'}(t) = \mathbf{A}\,\bm{u}(t) \qquad \bm{u}(0) = \bm{u}_0 The solution is :math:`\bm{u}(t) = \bm{u}_0 e^{\mathbf{A}\,t}`. Notice that we do not have a negative sign in front of the :math:`\bf{A}` matrix as we previously saw with the :math:`a` scalar ODE in :ref:`ode-stability`. But it is a decaying exponential that goes to zero when :math:`t` goes to infinity because the eigenvalues of :math:`\bf{A}` are negative. We get the starting equation for the numerical solution by applying :eq:`eq-implicitODE`. .. math:: \bm{u}_k = \bm{u}_{k+1} - h\,\mathbf{A}\,\bm{u}_{k+1} = (\mathbf{I} - h\,\mathbf{A})\,\bm{u}_{k+1} .. math:: \bm{u}_{k+1} = (\mathbf{I} - h\,\mathbf{A})^{-1}\,\bm{u}_k :label: eq-sysIBE :eq:`eq-sysIBE` is an implicit equation for solving systems of ODEs numerically, and is implemented in the ``IBEsysODE`` function listed below. Notice that the implicit approach requires solving an equation at every step. In this case, the equation is a linear system with a fixed step size, so the ``IBEsysODE`` function uses LU decomposition once and calls on |M|'s triangular solver with the left-divide operator to quickly find a solution in each loop iteration. If this were a non-linear problem, then numerical methods such as Newton's method or the secant method described in :ref:`fzero` would be needed, which could add significant computation. For this reason, implicit methods are not recommended for non-stiff or non-linear problems [SHEN15]_. :: function [t, u] = IBEsysODE(A, tspan, u0, n) % IBESYSODE - ODE solver using the Implicit Euler Backward % method for systems of ODEs. % % A - matrix defining the system coefficients % tspan - row vector [t0 tf] for initial and final time % y0 - initial value of y % n - number of y values to find h = (tspan(2) - tspan(1))/n; t = linspace(tspan(1), tspan(2), n); m = size(A, 1); u = zeros(m, n); u(:,1) = u0; [L, U, P] = lu(eye(m) - h*A); for k = 1:n-1 u(:,k+1) = U\(L\(P*u(:,k))); end :numref:`fig-IBEsys` shows the numerical solution that the ``IBEsysODE`` function found to our stiff system of ODEs. Since this is a fixed step size solution, the step size was decreased to 0.005 to accurately track the quickly changing :math:`e^{-49\,t}` terms in the solution. Of course, the small step is not needed once the :math:`e^{-49\,t}` term was sufficiently close to zero. Remember that there is a difference between accuracy and stability. A lack of stability causes the solution to be wrong as :math:`k \to \infty`. A lack of accuracy just means that the exact solution is not well modeled by the numerical solution, which is often solved by reducing the step size. :: >> A = [-25 24; 24 -25]; >> u0 = [1;2]; >> [t, u] = IBEsysODE(A, [0, 1], u0, 50); >> x = @(t) -0.5*exp(-49*t) + 1.5*exp(-t); >> y = @(t) 0.5*exp(-49*t) + 1.5*exp(-t); >> figure, hold on >> plot(t, u(1,:), 'o') >> plot(t, u(2,:), '+') >> plot(t, x(t), ':') >> plot(t, y(t), '--') >> hold off >> axis tight >> legend('u(1,:)', 'u(2,:)', 'x(t)', 'y(t)') >> xlabel('t') .. codeplot .. _fig-IBEsys: .. figure:: IBEsys.png :align: center :width: 80% With enough data points, the ``IBEsysODE`` function is able to to reasonably well track the rapidly changing solution of the :math:`e^{-49\,t}` terms. To evaluate the stability of the implicit backward Euler method, we start by using induction to express the steady state solution in terms of the initial value. .. math:: \bm{u}_k = \left[(\mathbf{I} - h\,\mathbf{A})^{-1}\right]^k\,u_0 The stability constraint is that :math:`\bm{u}_k \to 0` as :math:`k \to \infty`. So we need that .. math:: \norm{(I - h\,A)^{-1}}_2 < 1, \; \text{ for all }h. Since our :math:`\bf{A}` matrix is symmetric, we can use the property of :math:`l_2` matrix norms that :math:`\norm{\mathbf{A}}_2 = \max_i \abs{\lambda_i(\mathbf{A})}`, where :math:`\lambda_i(\mathbf{A})` means the :math:`i^{th}` eigenvalue of :math:`\bf{A}`. (See :ref:`matNorms` for more information of matrix norms.) .. math:: \max_i \abs{\lambda_i\left((I - h\,A)^{-1}\right)} < 1 Now we can make use of the eigenvalue property for symmetric matrices that :math:`\lambda_i(\mathbf{A}^{-1}) = \frac{1}{\lambda_i(\mathbf{A})}`. (See :ref:`eigProperties`) .. eq-eigInv .. math:: :label: eq-SysStable \begin{array}{rl} & \max_i \abs{\frac{1}{\lambda_i(I - h\,\mathbf{A})}} < 1 \\ \hfill & \min_i \abs{\lambda_i(I - h\,\mathbf{A})} = \min_i \abs{1 - h\,\lambda_i(\mathbf{A})} > 1, \end{array} :eq:`eq-SysStable` holds for any symmetric :math:`\bf{A}` matrix with negative eigenvalues. So the implicit backward Euler method is unconditionally stable for systems of ODEs. Unconditionally stable ODE solving algorithms remove the concern about having an unstable numerical solution because the step size is too large, but as we saw in :numref:`fig-IBEsys`, accuracy requirements may still require a reduced step size. The more complex and robust algorithms of |M|'s ODE solvers for stiff systems are able to adjust the step size enough to find accurate solutions without being excessively slowed by a stiff systems. MATLAB's ``ode15s`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |M| has several stiff ODE solvers. The ``ode15s`` should be tried first for stiff ODEs. The primary authors of |M|'s suite of ODE solvers, Shampine and Reichelt, describe ``ode15s`` as an implicit quasi-constant step size algorithm. They further give the following suggestion regarding picking which of their ODE solver to use [SHAMPINE97]_. :math:`\ldots` Except in special circumstances, ``ode45`` should be the code tried first. If there is reason to believe the problem to be stiff, or if the problem turns out to be unexpectedly difficult for ``ode45``, the ``ode15s`` code should be tried. :numref:`fig-ode15s_Sys` shows how ``ode15s`` handles the stiff system of ODEs used in the previous example. :: >> A = [-25 24; 24 -25]; >> fA = @(t, u) A*u; >> u0 = [1;2]; >> [t, y_ode15s] = ode15s(fA, [0 1], u0); >> t = t'; u = y_ode15s'; >> x = @(t) -0.5*exp(-49*t) + 1.5*exp(-t); >> y = @(t) 0.5*exp(-49*t) + 1.5*exp(-t); >> figure, hold on >> plot(t, u(1,:)', 'o') >> plot(t, u(2,:)', '+') >> plot(t, x(t), ':') >> plot(t, y(t), '--') >> hold off >> legend('u(1,:)', 'u(2,:)', 'x(t)', 'y(t)') >> xlabel('t'), ylabel('y') .. codeplot .. _fig-ode15s_Sys: .. figure:: ode15s_Sys.png :align: center :width: 80% The ``ode15s`` ODE solver is designed for stiff systems, so the less significant, but rapidly changing terms do not significantly reduce the step size. MATLAB's Suite of ODE Solvers -------------------------------- |M|'s ``ode45`` and ``ode15s`` are the nonstiff and stiff ODE solvers to use for most problems. If those do not give satisfactory results then ``ode23`` and ``ode113`` might be tried for nonstiff problems. For stiff problems, ``ode23s`` should be the next function tried after ``ode15s``. Additional documentation and research may be needed to select the best solver to use for some problems. The MathWorks' web pages have quite a bit of documentation on the suit of ODE solvers. The *MATLAB Guide* by Higham and Higham [HIGHAM17]_ has extended coverage |M|'s suite of ODE solvers. :numref:`tbl-ode-solvers` lists |M|'s suite of ODE solvers. .. _tbl-ode-solvers: .. table:: The |M| ODE solvers =========== ============== ========================================= Solver Problem type Description =========== ============== ========================================= ``ode45`` Nonstiff Explicit Runge-Kutta, orders 4 and 5 ``ode23`` Nonstiff Explicit Runge-Kutta, orders 2 and 3 ``ode113`` Nonstiff Explicit linear multistep, orders 1--13 ``ode15s`` Stiff Implicit linear multistep, orders 1--5 ``ode15i`` Fully implicit Implicit linear multistep, orders 1--5 ``ode23s`` Stiff Modified Rosenbrock pair, orders 2 and 3 ``ode23t`` Mildly stiff Implicit trapezoid rule, orders 2 and 3 ``ode23tb`` Stiff Implicit Runge-Kutta, orders 2 and 3 =========== ============== ========================================= .. topic:: Online Resources This `MathWorks web page `_ has some good tutorial videos about numerical solutions to differential equations. Please watch at least the first three videos (*Euler -- ODE1*, *Midpoint Method -- ODE2*, and *Classical Runge-Kutta -- ODE4*). Although, the ODE solvers in |M| are slightly more complex than the examples in the first three videos, they explain the basic concepts. |M| has several ODE solvers, they are reviewed on this `web page `_. .. raw:: latex \clearpage