{"id":990,"date":"2014-06-09T12:00:29","date_gmt":"2014-06-09T17:00:29","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=990"},"modified":"2017-02-24T20:11:50","modified_gmt":"2017-02-25T01:11:50","slug":"ordinary-differential-equations-stiffness","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2014\/06\/09\/ordinary-differential-equations-stiffness\/","title":{"rendered":"Ordinary Differential Equations, Stiffness"},"content":{"rendered":"<div class=\"content\">\r\n\r\n<!--introduction-->Stiffness is a subtle concept that plays an important role in assessing the effectiveness of numerical methods for ordinary differential equations. (This article is adapted from section 7.9, \"Stiffness\", in <i>Numerical Computing with MATLAB<\/i>.)\r\n\r\n<!--\/introduction-->\r\n<h3>Contents<\/h3>\r\n<div>\r\n<ul>\r\n\t<li><a href=\"#28b6c990-f835-4bcb-9b2f-e5cfdc365e17\">Stiffness<\/a><\/li>\r\n\t<li><a href=\"#1aac84a3-f82d-4514-a587-b944035e571e\">Flame example<\/a><\/li>\r\n\t<li><a href=\"#e3da3e55-d349-4278-a117-237604d1e31f\">Stiffness in action<\/a><\/li>\r\n\t<li><a href=\"#6319886b-f60f-407e-9f0b-25dc4add3f0a\">Stiff solver<\/a><\/li>\r\n\t<li><a href=\"#05f61f43-dec6-4460-a232-9a7d35f4a218\">Take a hike<\/a><\/li>\r\n\t<li><a href=\"#e3ba3f06-5176-4d6a-9676-030f7868d704\">Implicit methods<\/a><\/li>\r\n\t<li><a href=\"#c96a70b6-62fb-4dc3-a40f-b186664c6101\">Newton's method<\/a><\/li>\r\n\t<li><a href=\"#37f3fdb4-a386-44cd-adcd-adb6918fa2e2\">ode15s<\/a><\/li>\r\n\t<li><a href=\"#3542cccf-01b6-4ed1-990a-da088ff147e3\">ode23s<\/a><\/li>\r\n\t<li><a href=\"#ca2950f4-b98c-4c0a-8a2d-ed03ad291af4\">ode23t and ode23tb<\/a><\/li>\r\n\t<li><a href=\"#5a5763a8-74fc-461f-840a-cf4d411c2444\">Matrix computations<\/a><\/li>\r\n\t<li><a href=\"#2ac360e5-79e4-46e7-8644-abeb79b51b61\">References<\/a><\/li>\r\n\t<li><a href=\"#fd430004-cdd5-409f-8c7b-51031e4eb008\">Postscript<\/a><\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>Stiffness<a name=\"28b6c990-f835-4bcb-9b2f-e5cfdc365e17\"><\/a><\/h4>\r\nStiffness is a subtle, difficult, and important concept in the numerical solution of ordinary differential equations. It depends on the differential equation, the initial conditions, and the numerical method. Dictionary definitions of the word \"stiff\" involve terms like \"not easily bent,\" \"rigid,\" and \"stubborn.\" We are concerned with a computational version of these properties.\r\n<p style=\"margin-left: 3ex;\">A problem is stiff if the solution being sought varies slowly,\r\nbut there are nearby solutions that vary rapidly, so the numerical\r\nmethod must take small steps to obtain satisfactory results.<\/p>\r\nStiffness is an efficiency issue. If we weren't concerned with how much time a computation takes, we wouldn't be concerned about stiffness. Nonstiff methods can solve stiff problems; they just take a long time to do it.\r\n\r\nStiff ode solvers do more work per step, but take bigger steps. And we're not talking about modest differences here. For truly stiff problems, a stiff solver can be orders of magnitude more efficient, while still achieving a given accuracy.\r\n<h4>Flame example<a name=\"1aac84a3-f82d-4514-a587-b944035e571e\"><\/a><\/h4>\r\nA model of flame propagation provides an example. I learned about this example from Larry Shampine, one of the authors of the MATLAB ordinary differential equation suite. If you light a match, the ball of flame grows rapidly until it reaches a critical size. Then it remains at that size because the amount of oxygen being consumed by the combustion in the interior of the ball balances the amount available through the surface. The simple model is\r\n\r\n$$ \\dot{y} = y^2 - y^3 $$\r\n\r\n$$ y(0) = \\delta $$\r\n\r\n$$ 0 \\leq t \\leq {2 \\over \\delta} $$.\r\n\r\nThe scalar variable $y(t)$ represents the radius of the ball. The $y^2$ and $y^3$ terms come from the surface area and the volume. The critical parameter is the initial radius, $\\delta$, which is \"small.\" We seek the solution over a length of time that is inversely proportional to $\\delta$. You can see a dramatization by downloading <tt>flame.m,<\/tt>\u00a0<a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/37976-numerical-computing-with-matlab\">available here<\/a>.\r\n\r\nAt this point, I suggest that you start up MATLAB and actually run our examples. It is worthwhile to see them in action. I will start with <tt>ode45<\/tt>, the workhorse of the MATLAB ode suite. If $\\delta$ is not very small, the problem is not very stiff. Try $\\delta$ = 0.02 and request a relative error of $10^{-4}$.\r\n<pre class=\"codeinput\">   delta = 0.02;\r\n   F = @(t,y) y^2 - y^3;\r\n   opts = odeset(<span class=\"string\">'RelTol'<\/span>,1.e-4);\r\n   ode45(F,[0 2\/delta],delta,opts);\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/odes_stiff_01.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nWith no output arguments, <tt>ode45<\/tt> automatically plots the solution as it is computed. You should get a plot of a solution that starts at $y$ = 0.01, grows at a modestly increasing rate until $t$ approaches 50, which is $1\/\\delta$, then grows rapidly until it reaches a value close to 1, where it remains.\r\n<h4>Stiffness in action<a name=\"e3da3e55-d349-4278-a117-237604d1e31f\"><\/a><\/h4>\r\nNow let's see stiffness in action. Decrease $\\delta$ by more than two orders of magnitude. (If you run only one example, run this one.)\r\n<pre class=\"codeinput\">   delta = 0.0001;\r\n   ode45(F,[0 2\/delta],delta,opts);\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/odes_stiff_02.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nYou should see something like this plot, although it will take several seconds to complete the plot. If you get tired of watching the agonizing progress, click the stop button in the lower left corner of the window. Turn on zoom, and use the mouse to explore the solution near where it first approaches steady state. You should see something like the detail obtained with this <tt>axis<\/tt> command.\r\n<pre class=\"codeinput\">   axis([.995e4 1.03e4 0.9998 1.0002])\r\n   last_plot = getframe(gcf);\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/odes_stiff_03.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nNotice that <tt>ode45<\/tt> is doing its job. It's keeping the solution within $10^{-4}$ of its nearly constant steady state value. But it certainly has to work hard to do it. If you want an even more dramatic demonstration of stiffness, decrease $\\delta$ or decrease the tolerance to $10^{-5}$ or $10^{-6}$.\r\n\r\nThis problem is not stiff initially. It only becomes stiff as the solution approaches steady state. This is because the steady state solution is so \"rigid.\" Any solution near $y(t) = 1$ increases or decreases rapidly toward that solution. (I should point out that \"rapidly\" here is with respect to an unusually long time scale.)\r\n\r\nWhat can be done about stiff problems? You don't want to change the differential equation or the initial conditions, so you have to change the numerical method. Methods intended to solve stiff problems efficiently do more work per step, but can take much bigger steps. Stiff methods are <i>implicit<\/i>. At each step they use MATLAB matrix operations to solve a system of simultaneous linear equations that helps predict the evolution of the solution. For our flame example, the matrix is only 1 by 1, but even here, stiff methods do more work per step than nonstiff methods.\r\n<h4>Stiff solver<a name=\"6319886b-f60f-407e-9f0b-25dc4add3f0a\"><\/a><\/h4>\r\nLet's compute the solution to our flame example again, this time with <tt>ode23s<\/tt>. The \" <i>s<\/i> \" in the name is for \"stiff.\"\r\n<pre class=\"codeinput\">   ode23s(F,[0 2\/delta],delta,opts);\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/odes_stiff_04.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nHere is the zoom detail.\r\n<pre class=\"codeinput\">   axis([.995e4 1.03e4 0.9998 1.0002])\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/odes_stiff_05.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nYou can see that <i>ode23s<\/i> takes many fewer steps than <i>ode45<\/i>. This is actually an easy problem for a stiff solver. In fact, <i>ode23s<\/i> takes only 99 steps and uses just 412 function evaluations, while <i>ode45<\/i> takes 3040 steps and uses 20179 function evaluations. Stiffness even affects graphical output. The print files for the <i>ode45<\/i> figures are much larger than those for the <i>ode23s<\/i> figures.\r\n<h4>Take a hike<a name=\"05f61f43-dec6-4460-a232-9a7d35f4a218\"><\/a><\/h4>\r\nImagine you are returning from a hike in the mountains. You are in a narrow canyon with steep slopes on either side. An explicit algorithm would sample the local gradient to find the descent direction. But following the gradient on either side of the trail will send you bouncing back and forth across the canyon, as with #ode45#. You will eventually get home, but it will be long after dark before you arrive. An implicit algorithm would have you keep your eyes on the trail and anticipate where each step is taking you. It is well worth the extra concentration.\r\n<h4>Implicit methods<a name=\"e3ba3f06-5176-4d6a-9676-030f7868d704\"><\/a><\/h4>\r\nAll numerical methods for stiff odes are <i>implicit<\/i>. The simplest example is the <i>backward Euler<\/i> method, which involves finding the unknown $v$ in\r\n\r\n$$ v = y_n + h f(t_{n+1},v) $$\r\n\r\nand then setting $y_{n+1}$ equal to that $v$. This is usually a nonlinear system of equations whose solution requires at least an approximation to the Jacobian, the matrix of partial derivatives\r\n\r\n$$ J = {\\partial F \\over \\partial y} $$\r\n\r\nBy default the partial derivatives in the Jacobian are computed by finite differences. This can be quite costly in terms of function evaluations. If a procedure for computing the Jacobian is available, it can be provided. Or, if the sparsity pattern is known, it can be specified. The blocks in a Simulink diagram, for example, are only sparsely connected to each other. Specifying a sparse Jacobian initiates sparse linear equation solving.\r\n<h4>Newton's method<a name=\"c96a70b6-62fb-4dc3-a40f-b186664c6101\"><\/a><\/h4>\r\nNewton's method for computing the $v$ in the backward Euler method is an iteration. Start perhaps with $v^0 = y_n$. Then, for $k = 0,1,...$, solve the linear system of equations\r\n\r\n$$ (I - hJ) u^k = v^k - y_n - h f(t_{n+1},v^k) $$\r\n\r\nfor the correction $u^k$ . Set\r\n\r\n$$ v^{k+1} = v^k + u^k $$\r\n\r\nWhen you are satisfied that the $v^k$ have converged, let $y_{n+1}$ be the limit.\r\n\r\nStiff ODE solvers may not use Newton's method itself, but whatever method is used to find the solution, $y_{n+1}$, at the forward time step can ultimately be traced back to Newton's method.\r\n<h4>ode15s<a name=\"37f3fdb4-a386-44cd-adcd-adb6918fa2e2\"><\/a><\/h4>\r\n<tt>ode15s<\/tt> employs two variants of a method that is quite different from the single step methods that I've described so far in this series on ode solvers. Linear multistep methods save solution values from several time steps and use all of them to advance to the next step.\r\n\r\nActually, <tt>ode15s<\/tt> can be compared to the other multistep method in the suite, <tt>ode113<\/tt>. One saves values of the solution, $y_n$, while the other saves values of the function, $F(t_n,y_n)$. One includes the unknown value at the new time step $y_{n+1}$ in the formulation, thereby making it implicit, while the other does not. Both methods can vary the order as well as the step size. As their names indicate, <tt>ode15s<\/tt> allows the order to vary between 1 and 5, while <tt>ode113s<\/tt> allows the order to vary between 1 and 13.\r\n\r\nA property specified via <tt>odeset<\/tt> switches <tt>ode15s<\/tt> between two variants of a linear multistep method, BDF, Backward Differentiation Formulas, and NDF, Numerical Differentiation Formulas. BDFs were introduced for stiff odes in 1971 by Bill Gear. Gear's student, Linda Petzold, extended the ideas to DAEs, differential-algebraic equations, and produced DASSL, software whose successors are still in widespread use today. NDFs, which are the default for <tt>ode15s<\/tt>, include an additional term in the memory and consequently can take larger steps with the same accuracy, especially at lower order.\r\n<h4>ode23s<a name=\"3542cccf-01b6-4ed1-990a-da088ff147e3\"><\/a><\/h4>\r\n<tt>ode23s<\/tt> is a single step, implicit method described in the paper by Shampine and Reichelt referenced below. It uses a second order formula to advance the step and a third order formula to estimate the error. It recomputes the Jacobian with each step, thereby making it quite expensive in terms of function evaluations. If you can supply an analytic Jacobian then <tt>ode23s<\/tt> is a competitive choice.\r\n<h4>ode23t and ode23tb<a name=\"ca2950f4-b98c-4c0a-8a2d-ed03ad291af4\"><\/a><\/h4>\r\n<tt>ode23t<\/tt> and <tt>ode23tb<\/tt> are implicit methods based on the trapezoid rule and the second order BDF. The origins of the methods go back to the 1985 paper referenced below by a group at the old Bell Labs working on electronic device and circuit simulation. Mike Hosea and Larry Shampine made extensive modifications and improvements described in their 1996 paper when they implemented the methods in MATLAB.\r\n<h4>Matrix computations<a name=\"5a5763a8-74fc-461f-840a-cf4d411c2444\"><\/a><\/h4>\r\nStiff ODE solvers are not actually using MATLAB's iconic backslash operator on a full system of linear equations, but they are using its component parts, LU decomposition and solution of the resulting triangular systems.\r\n\r\nLet's look at the statistics generated by <tt>ode23<\/tt> when it solves the flame problem. We'll run it again, avoiding the plot by asking for output, but then ignoring the output, and just looking at the stats.\r\n<pre class=\"codeinput\">   opts = odeset(<span class=\"string\">'stats'<\/span>,<span class=\"string\">'on'<\/span>,<span class=\"string\">'reltol'<\/span>,1.e-4);\r\n   [~,~] = ode23s(F,[0 2\/delta],delta,opts);\r\n<\/pre>\r\n<pre class=\"codeoutput\">99 successful steps\r\n7 failed attempts\r\n412 function evaluations\r\n99 partial derivatives\r\n106 LU decompositions\r\n318 solutions of linear systems\r\n<\/pre>\r\nWe see that at every step <tt>ode23s<\/tt> is computing a Jacobian, finding the LU decomposition of a matrix involving that Jacobian, and then using L and U to solve three linear systems.\r\n\r\nNow how about the primary stiff solver, <tt>ode15s<\/tt>.\r\n<pre class=\"codeinput\">   [~,~] = ode15s(F,[0 2\/delta],delta,opts);\r\n<\/pre>\r\n<pre class=\"codeoutput\">140 successful steps\r\n39 failed attempts\r\n347 function evaluations\r\n2 partial derivatives\r\n53 LU decompositions\r\n342 solutions of linear systems\r\n<\/pre>\r\nWe see that <tt>ode15s<\/tt> takes more steps than <tt>ode23s<\/tt>, but requires only two Jacobians. It does only half as many LU decompositions, but then uses each LU decomposition for twice as many linear equation solutions.\r\n\r\nWe certainly can't draw any conclusions about the relative merits of these two solvers from this one example, especially since the Jacobian in this case is only a 1-by-1 matrix.\r\n<h4>References<a name=\"2ac360e5-79e4-46e7-8644-abeb79b51b61\"><\/a><\/h4>\r\nCleve Moler, <i>Numerical Computing with MATLAB<\/i>, Electronic Edition, MathWorks, <a href=\"https:\/\/www.mathworks.com\/moler\/index_ncm.html\">&lt;https:\/\/www.mathworks.com\/moler\/index_ncm.html<\/a>&gt;,\r\n\r\nPrint Edition, SIAM Revised Reprint, SIAM, 2008, 334 pp., <a href=\"https:\/\/bookstore.siam.org\/OT87\">&lt;https:\/\/bookstore.siam.org\/OT87<\/a>&gt;.\r\n\r\nLawrence F. Shampine and Mark W. Reichelt, \"The MATLAB ODE Suite\", SIAM Journal on Scientific Computing, 18 (1997), pp.1-22, <a href=\"https:\/\/www.mathworks.com\/help\/pdf_doc\/otherdocs\/ode_suite.pdf\">&lt;https:\/\/www.mathworks.com\/help\/pdf_doc\/otherdocs\/ode_suite.pdf<\/a>&gt;\r\n\r\nLawrence F. Shampine, Mark W. Reichelt, and Jacek A. Kierzenka, \"Solving Index-1 DAEs in MATLAB and Simulink\", SIAM Review, 41 (1999), pp. 538-552. <a href=\"http:\/\/epubs.siam.org\/doi\/abs\/10.1137\/S003614459933425X\">&lt;http:\/\/epubs.siam.org\/doi\/abs\/10.1137\/S003614459933425X<\/a>&gt;\r\n\r\nM. E. Hosea and L. F. Shampine, \"Analysis and Implementawtion of TR-BDF2\", Applied Numerical Mathematicss, 20 (1996), pp. 21-37, <a href=\"http:\/\/www.sciencedirect.com\/science\/article\/pii\/0168927495001158\">&lt;http:\/\/www.sciencedirect.com\/science\/article\/pii\/0168927495001158<\/a>&gt;\r\n\r\nR. E. Bank, W. M. Coughran, Jr., W. Fichtner., E. H. Grosse, D. J. Rose, and R. K. Smith, \"Transient simulation of silicon devices and circuits\", IEEE Transactions on Computer-Aided Design CAD-4 (1985), 4, pp. 436-451.\r\n<h4>Postscript<a name=\"fd430004-cdd5-409f-8c7b-51031e4eb008\"><\/a><\/h4>\r\nI want to repeat this plot because it represents the post on the <a href=\"https:\/\/blogs.mathworks.com\/cleve\">lead-in page at MATLAB Central<\/a>.\r\n<pre class=\"codeinput\">   imshow(last_plot.cdata,<span class=\"string\">'border'<\/span>,<span class=\"string\">'tight'<\/span>)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/odes_stiff_06.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/><script>\/\/ <![CDATA[\r\nfunction grabCode_206a3a4e7be54a3386edde19863cc2c2() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='206a3a4e7be54a3386edde19863cc2c2 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 206a3a4e7be54a3386edde19863cc2c2';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2014 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('\r\n\r\n\r\n\r\n<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\r\n\r\n\r\n\r\n\r\n\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }\r\n\/\/ ]]><\/script>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\"><a><span style=\"font-size: x-small; font-style: italic;\">Get\r\nthe MATLAB code<noscript>(requires JavaScript)<\/noscript><\/span><\/a><\/p>\r\nPublished with MATLAB\u00ae R2014a\r\n\r\n<\/div>\r\n<!--\r\n206a3a4e7be54a3386edde19863cc2c2 ##### SOURCE BEGIN #####\r\n%% Ordinary Differential Equations, Stiffness\r\n% Stiffness is a subtle concept that plays an important role in assessing\r\n% the effectiveness of numerical methods for ordinary differential equations.\r\n% (This article is adapted from section 7.9, \"Stiffness\", in _Numerical\r\n% Computing with MATLAB_.)\r\n%\r\n%% Stiffness\r\n% Stiffness is a subtle, difficult, and important concept in the numerical\r\n% solution of ordinary differential equations.  It depends on the\r\n% differential equation, the initial conditions, and the numerical method.\r\n% Dictionary definitions of the word \"stiff\" involve terms like\r\n% \"not easily bent,\" \"rigid,\" and \"stubborn.\"\r\n% We are concerned with a computational version of these properties.\r\n%\r\n% <html>\r\n%\r\n\r\n<p style=\"margin-left:3ex;\">\r\n% A problem is stiff if the solution being sought varies slowly,\r\n% but there are nearby solutions that vary rapidly, so the numerical\r\n% method must take small steps to obtain satisfactory results.\r\n%\r\n\r\n% <\/html>\r\n\r\n%%\r\n% Stiffness is an efficiency issue.  If we weren't concerned with how\r\n% much time a computation takes, we wouldn't be concerned about stiffness.\r\n% Nonstiff methods can solve stiff problems; they just take a long time\r\n% to do it.\r\n\r\n%%\r\n% Stiff ode solvers do more work per step, but take bigger steps.\r\n% And we're not talking about modest differences here.  For truly stiff\r\n% problems, a stiff solver can be orders of magnitude more efficient,\r\n% while still achieving a given accuracy.\r\n\r\n%% Flame example\r\n% A model of flame propagation provides an example.\r\n% I learned about this example from Larry Shampine, one of the authors\r\n% of the MATLAB ordinary differential equation suite.\r\n% If you light a match, the ball of flame grows\r\n% rapidly until it reaches a critical size.  Then it remains at that\r\n% size because the amount of oxygen being consumed by the combustion\r\n% in the interior of the ball balances the amount available through\r\n% the surface.  The simple model is\r\n%\r\n% $$ \\dot{y} = y^2 - y^3 $$\r\n%\r\n% $$ y(0) = \\delta $$\r\n%\r\n% $$ 0 \\leq t \\leq 2 \\leq \\delta $$.\r\n%\r\n% The scalar variable $y(t)$ represents the radius of the ball.  The\r\n% $y^2$ and $y^3$ terms come from the surface area and the volume.\r\n% The critical parameter is the initial radius, $\\delta$, which is \"small.\"\r\n% We seek the solution over a length of time that is inversely proportional\r\n% to $\\delta$.\r\n% You can see a dramatization by downloading |flame.m| from\r\n% <https:\/\/www.mathworks.com\/moler.htmlncmfilelist.html % https:\/\/www.mathworks.com\/moler.htmlncmfilelist.html>.\r\n\r\n%%\r\n% At this point, I suggest that you start up MATLAB and actually run our\r\n% examples.  It is worthwhile to see them in action.  I will start with\r\n% |ode45|, the workhorse of the MATLAB ode suite.\r\n% If $\\delta$ is not very small, the problem is not very stiff.\r\n% Try $\\delta$ = 0.02 and request a relative error of $10^{-4}$.\r\n%\r\ndelta = 0.02;\r\nF = @(t,y) y^2 - y^3;\r\nopts = odeset('RelTol',1.e-4);\r\node45(F,[0 2\/delta],delta,opts);\r\n\r\n%%\r\n% With no output arguments, |ode45| automatically plots the solution as it\r\n% is computed.  You should get a plot of a solution that starts at $y$ = 0.01,\r\n% grows at a modestly increasing rate until $t$ approaches 50, which is\r\n% $1\/\\delta$, then grows rapidly until it reaches a value close to 1,\r\n% where it remains.\r\n\r\n%% Stiffness in action\r\n% Now let's see stiffness in action.  Decrease $\\delta$ by more than two\r\n% orders of magnitude.  (If you run only one example, run this one.)\r\n\r\ndelta = 0.0001;\r\node45(F,[0 2\/delta],delta,opts);\r\n\r\n%%\r\n% You should see something like this plot, although it will take several\r\n% seconds to complete the plot.  If you get tired of watching the agonizing\r\n% progress, click the stop button in the lower left corner of the window.\r\n% Turn on zoom, and use the mouse to explore the solution near where it first\r\n% approaches steady state.  You should see something like the detail\r\n% obtained with this |axis| command.\r\n\r\naxis([.995e4 1.03e4 0.9998 1.0002])\r\nlast_plot = getframe(gcf);\r\n\r\n%%\r\n% Notice that |ode45| is doing its job.  It's keeping the solution within\r\n% $10^{-4}$ of its nearly constant steady state value.  But it certainly has\r\n% to work hard to do it.  If you want an even more dramatic demonstration of\r\n% stiffness, decrease $\\delta$ or decrease the tolerance to\r\n% $10^{-5}$ or $10^{-6}$.\r\n\r\n%%\r\n% This problem is not stiff initially.  It only becomes stiff as the solution\r\n% approaches steady state.  This is because the steady state solution is so\r\n% \"rigid.\"  Any solution near $y(t) = 1$ increases or decreases rapidly toward\r\n% that solution.  (I should point out that \"rapidly\" here is with respect\r\n% to an unusually long time scale.)\r\n%\r\n% What can be done about stiff problems?  You don't want to change the\r\n% differential equation or the initial conditions, so you have to change the\r\n% numerical method.  Methods intended to solve stiff problems efficiently do\r\n% more work per step, but can take much bigger steps.  Stiff methods are\r\n% _implicit_.  At each step they use MATLAB matrix operations to solve\r\n% a system of simultaneous linear equations that helps predict the evolution\r\n% of the solution.  For our flame example, the matrix is only 1 by 1, but\r\n% even here, stiff methods do more work per step than nonstiff methods.\r\n\r\n%% Stiff solver\r\n% Let's compute the solution to our flame example again, this time with\r\n% |ode23s|.  The \" _s_ \" in the name is for \"stiff.\"\r\n\r\node23s(F,[0 2\/delta],delta,opts);\r\n\r\n%%\r\n% Here is the zoom detail.\r\n\r\naxis([.995e4 1.03e4 0.9998 1.0002])\r\n\r\n%%\r\n% You can see that _ode23s_ takes many fewer steps than\r\n% _ode45_.  This is actually an easy problem for a stiff solver.  In fact,\r\n% _ode23s_ takes only 99 steps and uses just 412 function evaluations,\r\n% while _ode45_ takes 3040 steps and uses 20179 function evaluations.\r\n% Stiffness even affects graphical output.  The print files for the _ode45_\r\n% figures are much larger than those for the _ode23s_ figures.\r\n\r\n%% Take a hike\r\n% Imagine you are returning from a hike in the mountains.  You are in a narrow\r\n% canyon with steep slopes on either side.  An explicit algorithm would sample\r\n% the local gradient to find the descent direction.  But following the gradient\r\n% on either side of the trail will send you bouncing back and forth across\r\n% the canyon, as with #ode45#.  You will eventually get home, but it will be\r\n% long after dark before you arrive.  An implicit algorithm would have you\r\n% keep your eyes on the trail and anticipate where each step is taking you.\r\n% It is well worth the extra concentration.\r\n\r\n%% Implicit methods\r\n% All numerical methods for stiff odes are _implicit_.  The simplest example\r\n% is the _backward Euler_ method, which involves finding the unknown $v$ in\r\n%\r\n% $$ v = y_n + h f(t_{n+1},v) $$\r\n%\r\n% and then setting $y_{n+1}$ equal to that $v$.  This is usually a nonlinear\r\n% system of equations whose solution requires at least an approximation to\r\n% the Jacobian, the matrix of partial derivatives\r\n%\r\n% $$ J = {\\partial F \\over \\partial y} $$\r\n\r\n%%\r\n% By default the partial derivatives in the Jacobian are computed by finite\r\n% differences.  This can be quite costly in terms of function evaluations.\r\n% If a procedure for computing the Jacobian is available, it can be provided.\r\n% Or, if the sparsity pattern is known, it can be specified.  The blocks in\r\n% a Simulink diagram, for example, are only sparsely connected to each other.\r\n% Specifying a sparse Jacobian initiates sparse linear equation solving.\r\n\r\n%% Newton's method\r\n% Newton's method for computing the $v$ in the backward Euler method\r\n% is an iteration.  Start perhaps with $v^0 = y_n$.  Then,\r\n% for $k = 0,1,...$, solve the linear system of equations\r\n%\r\n% $$ (I - hJ) u^k = v^k - y_n - h f(t_{n+1},v^k) $$\r\n%\r\n% for the correction $u^k$ .  Set\r\n%\r\n% $$ v^{k+1} = v^k + u^k $$\r\n%\r\n% When you are satisfied that the $v^k$ have converged, let $y_{n+1}$ be\r\n% the limit.\r\n\r\n%%\r\n% Stiff ODE solvers may not use Newton's method itself, but whatever method is\r\n% used to find the solution, $y_{n+1}$, at the forward time step can ultimately\r\n% be traced back to Newton's method.\r\n\r\n%% ode15s\r\n% |ode15s| employs two variants of a method that is quite different from\r\n% the single step methods that I've described so far in this series on ode\r\n% solvers.  Linear multistep methods save solution values from several time\r\n% steps and use all of them to advance to the next step.\r\n\r\n%%\r\n% Actually, |ode15s| can be compared to the other multistep method in\r\n% the suite, |ode113|.  One saves values of the solution, $y_n$, while the\r\n% other saves values of the function, $F(t_n,y_n)$.  One includes the unknown\r\n% value at the new time step $y_{n+1}$ in the formulation, thereby making\r\n% it implicit, while the other does not.  Both methods can vary the order\r\n% as well as the step size.  As their names indicate, |ode15s| allows the\r\n% order to vary between 1 and 5, while |ode113s| allows the order to vary\r\n% between 1 and 13.\r\n\r\n%%\r\n% A property specified via |odeset| switches |ode15s| between two variants\r\n% of a linear multistep method, BDF, Backward Differentiation Formulas,\r\n% and NDF, Numerical Differentiation Formulas.\r\n% BDFs were introduced for stiff odes in 1971 by Bill Gear.  Gear's student,\r\n% Linda Petzold, extended the ideas to DAEs, differential-algebraic equations,\r\n% and produced DASSL, software whose successors are still in widespread use\r\n% today.\r\n% NDFs, which are the default for |ode15s|, include an additional term in\r\n% the memory and consequently can take larger steps with the same accuracy,\r\n% especially at lower order.\r\n\r\n%% ode23s\r\n% |ode23s| is a single step, implicit method described in the paper by\r\n% Shampine and Reichelt referenced below.  It uses a second order formula\r\n% to advance the step and a third order formula to estimate the error.\r\n% It recomputes the Jacobian with each step, thereby making it quite\r\n% expensive in terms of function evaluations.  If you can supply an\r\n% analytic Jacobian then |ode23s| is a competitive choice.\r\n\r\n%% ode23t and ode23tb\r\n% |ode23t| and |ode23tb| are implicit methods based on the trapezoid rule\r\n% and the second order BDF.  The origins of the methods go back to the 1985\r\n% paper referenced below by a group at the old Bell Labs working on electronic\r\n% device and circuit simulation.  Mike Hosea and Larry Shampine made extensive\r\n% modifications and improvements described in their 1996 paper when they\r\n% implemented the methods in MATLAB.\r\n\r\n%% Matrix computations\r\n% Stiff ODE solvers are not actually using MATLAB's iconic backslash\r\n% operator on a full system of linear equations, but they are using\r\n% its component parts, LU decomposition and solution of the resulting\r\n% triangular systems.\r\n\r\n%%\r\n% Let's look at the statistics generated by |ode23| when it solves\r\n% the flame problem.  We'll run it again, avoiding the plot by asking for\r\n% output, but then ignoring the output, and just looking at the stats.\r\n\r\nopts = odeset('stats','on','reltol',1.e-4);\r\n[~,~] = ode23s(F,[0 2\/delta],delta,opts);\r\n\r\n%%\r\n% We see that at every step |ode23s| is computing a Jacobian,\r\n% finding the LU decomposition of a matrix involving that Jacobian,\r\n% and then using L and U to solve three linear systems.\r\n\r\n%%\r\n% Now how about the primary stiff solver, |ode15s|.\r\n\r\n[~,~] = ode15s(F,[0 2\/delta],delta,opts);\r\n\r\n%%\r\n% We see that |ode15s| takes more steps than |ode23s|, but requires\r\n% only two Jacobians.  It does only half as many LU decompositions,\r\n% but then uses each LU decomposition for twice as many linear equation\r\n% solutions.\r\n\r\n%%\r\n% We certainly can't draw any conclusions about the relative merits of\r\n% these two solvers from this one example, especially since the Jacobian\r\n% in this case is only a 1-by-1 matrix.\r\n\r\n%% References\r\n% Cleve Moler, _Numerical Computing with MATLAB_, Electronic Edition,\r\n% MathWorks, <https:\/\/www.mathworks.com\/moler\/index_ncm.html % https:\/\/www.mathworks.com\/moler\/index_ncm.html>,\r\n%\r\n% Print Edition,\r\n% SIAM Revised Reprint, SIAM, 2008, 334 pp.,\r\n% <http:\/\/www.ec-securehost.com\/SIAM\/ot87.html % http:\/\/www.ec-securehost.com\/SIAM\/ot87.html>.\r\n\r\n%%\r\n% Lawrence F. Shampine and Mark W. Reichelt, \"The MATLAB ODE Suite\",\r\n% SIAM Journal on Scientific Computing, 18 (1997), pp.1-22,\r\n% <https:\/\/www.mathworks.com\/help\/pdf_doc\/otherdocs\/ode_suite.pdf % https:\/\/www.mathworks.com\/help\/pdf_doc\/otherdocs\/ode_suite.pdf>\r\n\r\n%%\r\n% Lawrence F. Shampine, Mark W. Reichelt, and Jacek A. Kierzenka,\r\n% \"Solving Index-1 DAEs in MATLAB and Simulink\", SIAM Review, 41 (1999),\r\n% pp. 538-552.  <http:\/\/epubs.siam.org\/doi\/abs\/10.1137\/S003614459933425X % http:\/\/epubs.siam.org\/doi\/abs\/10.1137\/S003614459933425X>\r\n\r\n%%\r\n% M. E. Hosea and L. F. Shampine, \"Analysis and Implementawtion of TR-BDF2\",\r\n% Applied Numerical Mathematicss, 20 (1996), pp. 21-37,\r\n% <http:\/\/www.sciencedirect.com\/science\/article\/pii\/0168927495001158 % http:\/\/www.sciencedirect.com\/science\/article\/pii\/0168927495001158>\r\n\r\n%%\r\n% R. E. Bank, W. M. Coughran, Jr., W. Fichtner., E. H. Grosse, D. J. Rose,\r\n% and R. K. Smith, \"Transient simulation of silicon devices and circuits\",\r\n% IEEE Transactions on Computer-Aided Design CAD-4 (1985), 4, pp. 436-451.\r\n\r\n%% Postscript\r\n% I want to repeat this plot because it represents the post on the\r\n% <https:\/\/blogs.mathworks.com\/cleve lead-in page at MATLAB Central>.\r\n\r\nimshow(last_plot.cdata,'border','tight')\r\n\r\n##### SOURCE END ##### 206a3a4e7be54a3386edde19863cc2c2\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/feature_image\/odes_stiff_06.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div><!--introduction-->Stiffness is a subtle concept that plays an important role in assessing the effectiveness of numerical methods for ordinary differential equations. (This article is adapted from section 7.9, \"Stiffness\", in <i>Numerical Computing with MATLAB<\/i>.)\r\n\r\n<!--\/introduction-->... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2014\/06\/09\/ordinary-differential-equations-stiffness\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":1028,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[24,16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/990"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=990"}],"version-history":[{"count":9,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/990\/revisions"}],"predecessor-version":[{"id":2357,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/990\/revisions\/2357"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media\/1028"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=990"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=990"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=990"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}