{"id":1220,"date":"2015-09-23T13:47:32","date_gmt":"2015-09-23T18:47:32","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=1220"},"modified":"2015-08-31T13:57:51","modified_gmt":"2015-08-31T18:57:51","slug":"ode-solver-selection-in-matlab","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2015\/09\/23\/ode-solver-selection-in-matlab\/","title":{"rendered":"ODE Solver Selection in MATLAB"},"content":{"rendered":"\r\n<div class=\"content\"><!--introduction--><p>Today, I'd like to welcome Josh Meyer as this week's guest blogger. Josh works on the Documentation team here at MathWorks, where he writes and maintains some of the MATLAB Mathematics documentation. In this post, Josh provides a bit of advice on how to choose which ODE solver to use. Over to you, Josh...<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#328b2ed1-8935-4748-a468-062c22e23d8d\">Initial Value Problems<\/a><\/li><li><a href=\"#f7226a5d-a308-4a16-9b07-319010cf4c56\">Example: Euler's Method<\/a><\/li><li><a href=\"#06c88ca5-873d-40e7-95de-a2fa6e6dfc9f\">Improving on Euler's Method<\/a><\/li><li><a href=\"#25a5e77e-b415-4efd-8aeb-dd89e837f990\">Stiff Differential Equations<\/a><\/li><li><a href=\"#71d35ade-c8c2-4b37-8b6c-1bd9cfc77737\">Solver Recommendations<\/a><\/li><li><a href=\"#229371be-9d41-4619-ba29-09b81b57beac\">Example 1: Damped Pendulum<\/a><\/li><li><a href=\"#07f7ffa9-f2ba-477c-ab19-6f6c2d20e4e1\">Example 2: van der Pol Oscillator<\/a><\/li><li><a href=\"#15751569-d461-4f8d-afa6-7cd24363a220\">Summary<\/a><\/li><li><a href=\"#afce60fc-d932-421c-b7f2-48491eba1b82\">Comments<\/a><\/li><li><a href=\"#e197fc18-9fc4-4893-8c75-1316a54d1467\">Further Reading<\/a><\/li><\/ul><\/div><h4>Initial Value Problems<a name=\"328b2ed1-8935-4748-a468-062c22e23d8d\"><\/a><\/h4><p>There are 7 ordinary differential equation initial value problem solvers in MATLAB:<\/p><div><ul><li><tt>ode45<\/tt><\/li><li><tt>ode23<\/tt><\/li><li><tt>ode113<\/tt><\/li><li><tt>ode15s<\/tt><\/li><li><tt>ode23s<\/tt><\/li><li><tt>ode23t<\/tt><\/li><li><tt>ode23tb<\/tt><\/li><\/ul><\/div><p><i>(note that <tt>ode15i<\/tt> is left out of this discussion because it solves its own class of initial value problems: fully implicit ODEs of the form $f(t,y,y') = 0$)<\/i><\/p><p>To choose between the solvers, it's first necessary to understand why one solver might be better than another for a given problem.<\/p><p>The ODE solvers in MATLAB all work on initial value problems of the form,<\/p><p>$$y' = f \\left( t,y \\right)$$<\/p><p>where $y' = dy\/dt$. There is also a more general form,<\/p><p>$$ M(t,y) y' = f \\left( t,y \\right)$$<\/p><p>where $M(t,y)$ is referred to as the <i>mass matrix<\/i>.<\/p><p>Starting with the initial conditions $y_0$, and a period of time over which the answer is to be obtained $(t_0,t_f)$, the solution is obtained iteratively by using the results of previous steps according to the solver's algorithm. At the first such step, the initial conditions provide the necessary information that allows the integration to proceed. The final result is that the ODE solver returns a vector of time steps $t_0,t_1,...,t_f$ as well as the corresponding solution at each time step $y_0,y_1,...,y_f$.<\/p><p>Theoretically, this numerical solution technique is possible because of the connection between differential equations and integrals provided by the fundamental theorem of calculus:<\/p><p>$$y(t + h) = y(t) + \\int_t^{t+h} f\\left( s,y(s) \\right) ds$$<\/p><p>The problem of calculating $y(t+h)$ becomes a question of how to approximate the integral on the right hand side. This is where different solvers come in. Each different solver evaluates the integral using different numerical techniques, and each solver makes trade-offs between efficiency and accuracy.<\/p><h4>Example: Euler's Method<a name=\"f7226a5d-a308-4a16-9b07-319010cf4c56\"><\/a><\/h4><p>Euler's method is a simple ODE solver, but it provides an illustration of the trade-offs between efficiency and accuracy in an ODE solver algorithm. Suppose you want to solve<\/p><p>$$y' = f(t,y) = 2t$$<\/p><p>over the time span $[0,3]$ using the initial condition $y_0 = 0$. Each step of Euler's method is computed with<\/p><p>$$\\begin{array}{cl} y_{n+1} &amp;= y_n + h f \\left(t_n,y_n\\right)\\\\ t_{n+1}&amp;=\r\nt_n + h\\end{array}$$<\/p><p>Using $h=1$, the solution requires just three steps:<\/p><p>$$\\begin{array}{cl} y_1 &amp;= y_0 + f\\left(t_0,y_0\\right)=0\\\\y_2 &amp;= y_1 +\r\nf\\left(t_1,y_1 \\right)=2\\\\y_3 &amp;= y_2 +  f\\left(t_2,y_2 \\right)=6\r\n\\end{array}$$<\/p><p>... But is it accurate?<\/p><p>Not really. The exact solution to this equation is<\/p><p>$$y(t) = t^2$$<\/p><p>Reducing the step size $h$ can improve the accuracy of the answer a bit, but it also requires more steps to achieve the solution. To see this, the below code solves this problem using Euler's method and compares the answer to the analytic solution for several different values of $h$.<\/p><pre class=\"codeinput\">clear, clc\r\nh = 1;\r\ntspan = [0 3];\r\nf = @(t,y) 2*t;\r\ndydt(1) = 0;\r\nt(1) = 0;\r\n\r\ny = @(t) t.^2;\r\nx = linspace(0,tspan(end));\r\nfigure\r\nplot(x,y(x))\r\nxlabel(<span class=\"string\">'t'<\/span>), ylabel(<span class=\"string\">'y(t)'<\/span>)\r\nhold <span class=\"string\">on<\/span>\r\n\r\n<span class=\"keyword\">while<\/span> h &gt; 0.1\r\n    <span class=\"keyword\">for<\/span> k = 2:tspan(end)\/h+1\r\n        dydt(k) = dydt(k-1) + h*f(t(k-1),dydt(k-1));\r\n        t(k) = t(k-1) + h;\r\n    <span class=\"keyword\">end<\/span>\r\n\r\n    plot(t,dydt,<span class=\"string\">'-o'<\/span>)\r\n    h = h\/2;\r\n<span class=\"keyword\">end<\/span>\r\nlegend(<span class=\"string\">'Exact Solution'<\/span>, <span class=\"string\">'h=1'<\/span>,<span class=\"string\">'h=0.5'<\/span>,<span class=\"string\">'h=0.25'<\/span>,<span class=\"string\">'h=0.125'<\/span>,<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'Location'<\/span>,<span class=\"string\">'NorthWest'<\/span>)\r\ntitle(<span class=\"string\">'Solution of y''=2t using Euler''s method with several step sizes'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/OrdinaryDifferentialEquationsInMATLAB_01.png\" alt=\"\"> <h4>Improving on Euler's Method<a name=\"06c88ca5-873d-40e7-95de-a2fa6e6dfc9f\"><\/a><\/h4><p>Using smaller and smaller step sizes turns out to not be a good idea, since the algorithm loses efficiency. For any reasonable problem such a solver would be very slow. Also, Euler's method has a few inherent problems. Since the slope of $y$ is evaluated only once at the beginning of each interval, this solver only produces exact answers for constant functions. There is also no way to estimate the error, so the solver needs to use fixed step sizes.<\/p><p>So, one way to improve on Euler's method is to evaluate $y'$ more often in each step. This provides intermediate slopes that give a better idea of what the function is doing within each interval, allowing the solver to produce exact answers for higher order problems. For example, if you add an evaluation of the slope halfway across each interval to Euler's method, then the result is called the <i>midpoint rule<\/i>, which produces exact integrations for linear functions:<\/p><p>$$\\begin{array}{cl} s_1 &amp;= f(t_n,y_n)\\\\ s_2 &amp;= f\\left( t_n + \\frac{h}{2},\r\ny_n + \\frac{h}{2}s_1\\right)\\\\y_{n+1} &amp;= y_n + hs_2\\\\t_{n+1} &amp;=\r\nt_n+h\\end{array}$$<\/p><p>If you evaluate the slope four times in each interval, you get the <i>classical Runge-Kutta<\/i> algorithm (a.k.a. <i>RK4<\/i>), which is a piece of the <tt>ode45<\/tt> algorithm. This algorithm produces exact integrations for cubic functions (and if $f$ is only a function of $t$, then $s_2=s_3$ and this is the same as <i>Simpson's rule<\/i> for quadrature):<\/p><p>$$\\begin{array}{cl} s_1 &amp;= f(t_n,y_n)\\\\s_2 &amp;=\r\nf\\left(t_n+\\frac{h}{2},y_n+\\frac{h}{2}s_1\\right)\\\\s_3 &amp;=\r\nf\\left(t_n+\\frac{h}{2},y_n+\\frac{h}{2}s_2\\right)\\\\ s_4 &amp;=\r\nf\\left(t_n+h,y_n+hs_3\\right)\\\\y_{n+1} &amp;= y_n +\r\n\\frac{h}{6}\\left(s_1+2s_2+2s_3+s_4\\right)\\\\t_{n+1} &amp;= t_n +\r\nh\\end{array}$$<\/p><p>Runge-Kutta algorithms are all <i>single-step<\/i> solvers, since each step only depends on the result of the previous step. <tt>ode45<\/tt>, <tt>ode23<\/tt>, <tt>ode23s<\/tt>, <tt>ode23t<\/tt>, and <tt>ode23tb<\/tt> all employ single-step algorithms. Multi-step algorithms, such as those employed by <tt>ode113<\/tt> and <tt>ode15s<\/tt>, use the results of several past steps.<\/p><p>Sophisticated ODE solvers, like the ones in MATLAB, also estimate the error in each step to determine how big the next step size should be. This is another improvement over the fixed step sizes used above, since a solver that does more work per step is able to compensate by taking steps of varying size. The error estimate used to determine the step size is typically obtained by comparing the results of two different methods. MATLAB's ODE solvers follow a naming convention that reveals information about which methods they use. <tt>ode45<\/tt> compares the results of a 4th-order Runge-Kutta method and a 5th-order Runge-Kutta method to determine the error. Similarly, <tt>ode23<\/tt> uses a 2nd-order and 3rd-order Runge-Kutta comparison. So, in general, the smaller the number <tt>odeNN<\/tt>, the looser the solver's error tolerance is.<\/p><p>It should be no surprise, then, that <tt>ode45<\/tt> obtains a very accurate answer for the equation we solved before with Euler's method. <tt>ode45<\/tt> is MATLAB's general purpose ODE solver, and it is the first solver you should use for most problems.<\/p><pre class=\"codeinput\">y = @(t) t.^2;\r\nx = linspace(0,3);\r\nfigure\r\nplot(x,y(x))\r\nxlabel(<span class=\"string\">'t'<\/span>), ylabel(<span class=\"string\">'y(t)'<\/span>)\r\nhold <span class=\"string\">on<\/span>\r\n[t,y] = ode45(@(t,y) 2*t, [0 3], 0);\r\nplot(t,y,<span class=\"string\">'o'<\/span>)\r\nxlabel(<span class=\"string\">'t'<\/span>), ylabel(<span class=\"string\">'y(t)'<\/span>)\r\ntitle(<span class=\"string\">'Solution of y''=2t using ode45'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/OrdinaryDifferentialEquationsInMATLAB_02.png\" alt=\"\"> <h4>Stiff Differential Equations<a name=\"25a5e77e-b415-4efd-8aeb-dd89e837f990\"><\/a><\/h4><p>For some ODE problems, the step size taken by the solver is forced down to an unreasonably small level in comparison to the interval of integration, even in a region where the solution curve is smooth. These step sizes can be so small that traversing a short time interval might require millions of evaluations. This can lead to the solver failing the integration, but even if it succeeds it will take a very long time to do so.<\/p><p>Equations that cause this behavior in ODE solvers are said to be <i>stiff<\/i>. This is a nod to the fact that the equations are stubborn and not easily evaluated with numerical techniques. The problem that stiff ODEs pose is that explicit solvers (such as <tt>ode45<\/tt>) are untenably slow in achieving a solution. This is why <tt>ode45<\/tt> is classified as a nonstiff solver along with <tt>ode23<\/tt> and <tt>ode113<\/tt>. These solvers all struggle to integrate stiff equations.<\/p><p>Equation stiffness resists a precise definition, because there are several factors that cause it. Stiffness results from a combination of the specific equations, the ODE solver being used, the initial conditions, and the error tolerance used by the solver. The following statements about stiffness, attributed to Lambert [6], are exhibited by many examples of stiff ODEs, but counterexamples also exist, so they are not true definitions of stiffness:<\/p><div><ol><li>A linear constant coefficient system is stiff if all of its eigenvalues have negative real part and the stiffness ratio [of the largest and smallest eigenvalues] is large.<\/li><li>Stiffness occurs when the mathematical problem is stable, and yet stability requirements, rather than those of accuracy, severely constrain the step length.<\/li><li>Stiffness occurs when some components of the solution decay much more rapidly than others.<\/li><\/ol><\/div><p>A common theme among these statements is that stiffness can result from a difference in scaling somewhere in the problem. This difference in scale (for example, if the Jacobian $J = \\partial f_n\/\\partial y_i$ has a large ratio of negative eigenvalues) constrains the step size that the solver can take in performing the integration. Tiny step sizes become necessary in order to preserve any notion of error tolerance or stability in the solution.<\/p><p>For example, equations describing chemical reactions frequently display stiffness, since it is common for components of the solution to vary on drastically different time scales (reactions occurring at the same time that are both very slow and very fast).<\/p><p>However, there are solvers specifically designed to work on stiff ODEs. Solvers that are designed for stiff problems typically do more work per step, and the pay-off is that they are able to take much larger steps and enjoy improved numerical stability compared to the nonstiff solvers. Stiff solvers are implicit, because the computation of $y$ requires the use of linear algebra to solve systems of linear equations. The Jacobian is used to estimate the local behavior of the ODE as the integration proceeds, so supplying the analytical Jacobian can improve the performance of MATLAB's stiff ODE solvers.<\/p><p>This is just a cursory treatment of stiffness, because it is a complex topic. See <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2014\/06\/09\/ordinary-differential-equations-stiffness\/\">Ordinary Differential Equations: Stiffness<\/a> for a more in-depth look.<\/p><p>To summarize, the nonstiff solvers in MATLAB are:<\/p><div><ul><li><tt>ode45<\/tt><\/li><li><tt>ode23<\/tt><\/li><li><tt>ode113<\/tt><\/li><\/ul><\/div><p>The stiff solvers are (when <tt>ode45<\/tt> is slow):<\/p><div><ul><li><tt>ode15s<\/tt><\/li><li><tt>ode23s<\/tt><\/li><li><tt>ode23t<\/tt><\/li><li><tt>ode23tb<\/tt><\/li><\/ul><\/div><p>It should be noted that nonstiff solvers do <i>work<\/i> on stiff problems, it is just that they are exceptionally slow. Similarly, solvers designed for stiff problems can work on nonstiff problems, but since they do more work per step they are less efficient than their nonstiff counterparts when that extra work isn't necessary. So equation stiffness is a matter of solver efficiency, and the goal is to strike the right balance between accuracy of the solution and work done in each step by the solver.<\/p><h4>Solver Recommendations<a name=\"71d35ade-c8c2-4b37-8b6c-1bd9cfc77737\"><\/a><\/h4><p>The following recommendations are adapted from the <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/math\/choose-an-ode-solver.html\">MATLAB Mathematics documentation<\/a> :<\/p><div><ul><li><tt>ode45<\/tt> is MATLAB's general purpose single-step ODE solver. This should   be the first solver you use for most problems.<\/li><\/ul><\/div><p>For <b>nonstiff<\/b> problems:<\/p><div><ul><li><tt>ode23<\/tt> is another single-step solver that can be more efficient than   <tt>ode45<\/tt> if the problem permits a crude error tolerance. This looser   error tolerance can also accommodate some mildly stiff problems.<\/li><li><tt>ode113<\/tt> is a multi-step solver, and is preferred over <tt>ode45<\/tt> if the   function is expensive to evaluate, or for smooth problems where high   precision is required. For example, <tt>ode113<\/tt> excels with orbital   dynamics and celestial mechanics problems.<\/li><\/ul><\/div><p>For <b>stiff<\/b> problems (where <tt>ode45<\/tt> is slow):<\/p><div><ul><li><tt>ode15s<\/tt> is a multi-step solver that is MATLAB's general purpose solver   for stiff problems. Use <tt>ode15s<\/tt> if <tt>ode45<\/tt> fails or struggles to   complete the integration in a reasonable amount of time. <tt>ode15s<\/tt> is   also the primary solver for DAEs, which are identified as ODEs with a   singular mass matrix.<\/li><li>For stiff problems with crude error tolerances, <tt>ode23s<\/tt>, <tt>ode23t<\/tt>, and   <tt>ode23tb<\/tt> provide more efficient alternatives to <tt>ode15s<\/tt> since they   are single-step solvers. The efficiency of <tt>ode23s<\/tt> can be   significantly improved by providing the Jacobian, since <tt>ode23s<\/tt>   evaluates the Jacobian in each step.<\/li><li><tt>ode23s<\/tt> only works on ODEs with a mass matrix if the mass matrix is   constant (not time- or state-dependent).<\/li><li><tt>ode15s<\/tt> and <tt>ode23t<\/tt> are the only solvers that solve DAEs of index 1.<\/li><\/ul><\/div><p>Here is a graphic that captures the basic recommendations. In most cases, the only choice in solver you will need to make is to use <tt>ode15s<\/tt> instead of <tt>ode45<\/tt>.<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/selector.png\" alt=\"\"> <\/p><h4>Example 1: Damped Pendulum<a name=\"229371be-9d41-4619-ba29-09b81b57beac\"><\/a><\/h4><p>The equation of motion for a damped pendulum is,<\/p><p>$$\\ddot{\\theta} = -\\frac{b}{m}\\dot{\\theta}-\\frac{mg}{L\\left(m-2b\\right)}\r\n\\sin \\theta$$<\/p><p>where $g$ is the gravitational constant, $m$ the mass of the bob, $L$ the length of the string, and $b$ is a damping coefficient. The goal is to solve for $\\theta$, the angle that the pendulum deviates from the vertical, and $\\theta '$, the rate at which the angle changes.<\/p><p>Some natural initial conditions would be $\\theta_0 = \\pi\/4$ and $\\theta '_0 = 0$, indicating that you lift the pendulum up to a 45 degree angle before letting go, and it has no initial angular velocity. Due to the damping coefficient, you would expect the pendulum to slowly lose momentum and go back down to rest.<\/p><p>The file <tt>pendulumODE.m<\/tt> reformulates the problem as a coupled system of first-order ODEs:<\/p><p>$$\\begin{array}{cl} y_1' &amp;= y_2\\\\y_2' &amp;= -\\frac{b}{m} y_2\r\n-\\frac{mg}{L(m-2b)}sin(y_1)\\end{array}$$<\/p><p>then solves using <tt>ode45<\/tt>, <tt>ode15s<\/tt>, <tt>ode23<\/tt>, and <tt>ode113<\/tt>. The solutions for $y_1 = \\theta$ are plotted, and the file returns the stats for each solver. As is always the case when displaying execution times, \"the timings displayed can vary\".<\/p><pre class=\"language-matlab\">\r\n<span class=\"keyword\">function<\/span> pendulumODE\r\nopts = odeset(<span class=\"string\">'stats'<\/span>,<span class=\"string\">'on'<\/span>);\r\ntspan = [0 25];\r\ny0 = [pi\/4, 0];\r\n\r\ndisp(<span class=\"string\">' '<\/span>), disp(<span class=\"string\">'Stats for ode45:'<\/span>)\r\ntic, [t1,y1] = ode45(@pendode, tspan, y0, opts);\r\ntoc\r\n\r\ndisp(<span class=\"string\">' '<\/span>), disp(<span class=\"string\">'Stats for ode15s:'<\/span>)\r\ntic, [t2,y2] = ode15s(@pendode, tspan, y0, opts);\r\ntoc\r\n\r\ndisp(<span class=\"string\">' '<\/span>), disp(<span class=\"string\">'Stats for ode23:'<\/span>)\r\ntic, [t3,y3] = ode23(@pendode, tspan, y0, opts);\r\ntoc\r\n\r\ndisp(<span class=\"string\">' '<\/span>), disp(<span class=\"string\">'Stats for ode113:'<\/span>)\r\ntic, [t4,y4] = ode113(@pendode, tspan, y0, opts);\r\ntoc\r\n\r\nfigure\r\nsubplot(2,2,1), plot(t1,y1(:,1),<span class=\"string\">'-o'<\/span>), xlim([0 25]), title(<span class=\"string\">'ode45'<\/span>)\r\nsubplot(2,2,2), plot(t2,y2(:,1),<span class=\"string\">'-o'<\/span>), xlim([0 25]), title(<span class=\"string\">'ode15s'<\/span>)\r\nsubplot(2,2,3), plot(t3,y3(:,1),<span class=\"string\">'-o'<\/span>), xlim([0 25]), title(<span class=\"string\">'ode23'<\/span>)\r\nsubplot(2,2,4), plot(t4,y4(:,1),<span class=\"string\">'-o'<\/span>), xlim([0 25]), title(<span class=\"string\">'ode113'<\/span>)\r\n\r\n    <span class=\"keyword\">function<\/span> dy2dt2 = pendode(t,y)\r\n        g = 9.8; <span class=\"comment\">%m\/s^2<\/span>\r\n        m = 1;   <span class=\"comment\">% Mass of bob<\/span>\r\n        L = 2;   <span class=\"comment\">% Length of pendulum in meters<\/span>\r\n        b = 0.2; <span class=\"comment\">% Damping coefficient<\/span>\r\n        dy2dt2 = [y(2); -b\/m*y(2)-g\/L*sin(y(1))];\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n\r\n\r\n<\/pre><pre class=\"codeinput\">pendulumODE\r\n<\/pre><pre class=\"codeoutput\"> \r\nStats for ode45:\r\n75 successful steps\r\n0 failed attempts\r\n451 function evaluations\r\nElapsed time is 0.011970 seconds.\r\n \r\nStats for ode15s:\r\n183 successful steps\r\n9 failed attempts\r\n315 function evaluations\r\n1 partial derivatives\r\n31 LU decompositions\r\n311 solutions of linear systems\r\nElapsed time is 0.081467 seconds.\r\n \r\nStats for ode23:\r\n263 successful steps\r\n36 failed attempts\r\n898 function evaluations\r\nElapsed time is 0.015948 seconds.\r\n \r\nStats for ode113:\r\n159 successful steps\r\n3 failed attempts\r\n322 function evaluations\r\nElapsed time is 0.018056 seconds.\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/OrdinaryDifferentialEquationsInMATLAB_03.png\" alt=\"\"> <p>The solvers all perform well, but the damped pendulum is a good example of a nonstiff problem where <tt>ode45<\/tt> performs nicely. In this case <tt>ode15s<\/tt> needs to do extra work in order to achieve an inferior solution.<\/p><h4>Example 2: van der Pol Oscillator<a name=\"07f7ffa9-f2ba-477c-ab19-6f6c2d20e4e1\"><\/a><\/h4><p>The van der Pol Oscillator equation becomes stiff in certain intervals when the nonlinear parameter $\\mu$ is large:<\/p><p>$$\\ddot{y} - \\mu \\left(1-y^2\\right)\\dot{y}+y=0$$<\/p><p>The nonlinearity of this equation is contained entirely in the term that involves $\\mu$: notice that if $\\mu=0$, the equation reduces to that of a simple harmonic oscillator, which has regular periodic behavior.<\/p><p>Attempting to solve this equation using <tt>ode45<\/tt> is met with severe resistance, requiring millions of evaluations and 30+ minutes of execution (I stopped execution after 35 minutes). Since the problem is clearly stiff, this example compares the stiff solvers.<\/p><p>The file <tt>vanderpolODE.m<\/tt> finds the solution for $\\mu=1000$ using <tt>ode15s<\/tt>, <tt>ode23s<\/tt>, <tt>ode23t<\/tt>, and <tt>ode23tb<\/tt>. The function file <tt>vdp1000.m<\/tt> ships with MATLAB and encodes this equation as a coupled system of first-order ODEs:<\/p><p>$$\\begin{array}{cl}y_1' &amp;= y_2\\\\y_2' &amp;= \\mu \\left( 1-y_1^2 \\right)y_2 -\r\ny_1\\end{array}$$<\/p><p>The Jacobian is supplied to assist the solvers, and its use is reflected in the number of partial derivative evaluations.<\/p><pre class=\"language-matlab\">\r\n<span class=\"keyword\">function<\/span> vanderpolODE\r\nopts = odeset(<span class=\"string\">'stats'<\/span>,<span class=\"string\">'on'<\/span>,<span class=\"string\">'Jacobian'<\/span>,@J);\r\ntspan = [0 3000];\r\ny0 = [2 0];\r\n\r\ndisp(<span class=\"string\">' '<\/span>), disp(<span class=\"string\">'Stats for ode15s:'<\/span>)\r\ntic, [t1,y1] = ode15s(@vdp1000, tspan, y0, opts);\r\ntoc\r\n\r\ndisp(<span class=\"string\">' '<\/span>), disp(<span class=\"string\">'Stats for ode23s:'<\/span>)\r\ntic, [t2,y2] = ode23s(@vdp1000, tspan, y0, opts);\r\ntoc\r\n\r\ndisp(<span class=\"string\">' '<\/span>), disp(<span class=\"string\">'Stats for ode23t:'<\/span>)\r\ntic, [t3,y3] = ode23t(@vdp1000, tspan, y0, opts);\r\ntoc\r\n\r\ndisp(<span class=\"string\">' '<\/span>), disp(<span class=\"string\">'Stats for ode23tb:'<\/span>)\r\ntic, [t4,y4] = ode23tb(@vdp1000, tspan, y0, opts);\r\ntoc\r\n\r\nfigure\r\nsubplot(2,2,1), plot(t1,y1(:,1),<span class=\"string\">'-o'<\/span>), ylim([-4 4]), title(<span class=\"string\">'ode15s'<\/span>)\r\nsubplot(2,2,2), plot(t2,y2(:,1),<span class=\"string\">'-o'<\/span>), title(<span class=\"string\">'ode23s'<\/span>)\r\nsubplot(2,2,3), plot(t3,y3(:,1),<span class=\"string\">'-o'<\/span>), title(<span class=\"string\">'ode23t'<\/span>)\r\nsubplot(2,2,4), plot(t4,y4(:,1),<span class=\"string\">'-o'<\/span>), title(<span class=\"string\">'ode23tb'<\/span>)\r\n\r\n    <span class=\"keyword\">function<\/span> dfdy = J(t,y)\r\n        MU = 1000; <span class=\"comment\">% Nonlinear coefficient<\/span>\r\n        dfdy = [         0                  1\r\n            -2*MU*y(1)*y(2)-1    MU*(1-y(1)^2) ];\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n\r\n<\/pre><pre class=\"codeinput\">vanderpolODE\r\n<\/pre><pre class=\"codeoutput\"> \r\nStats for ode15s:\r\n591 successful steps\r\n225 failed attempts\r\n1749 function evaluations\r\n45 partial derivatives\r\n289 LU decompositions\r\n1747 solutions of linear systems\r\nElapsed time is 0.192955 seconds.\r\n \r\nStats for ode23s:\r\n741 successful steps\r\n13 failed attempts\r\n2251 function evaluations\r\n741 partial derivatives\r\n754 LU decompositions\r\n2262 solutions of linear systems\r\nElapsed time is 0.092824 seconds.\r\n \r\nStats for ode23t:\r\n776 successful steps\r\n94 failed attempts\r\n2014 function evaluations\r\n36 partial derivatives\r\n294 LU decompositions\r\n2012 solutions of linear systems\r\nElapsed time is 0.218996 seconds.\r\n \r\nStats for ode23tb:\r\n573 successful steps\r\n93 failed attempts\r\n2816 function evaluations\r\n44 partial derivatives\r\n269 LU decompositions\r\n3415 solutions of linear systems\r\nElapsed time is 0.202237 seconds.\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/OrdinaryDifferentialEquationsInMATLAB_04.png\" alt=\"\"> <p>The plots are of the solutions for $y_1$. For this problem, <tt>ode23s<\/tt> executes quickest and with the least number of failed steps. The supplied Jacobian greatly assists <tt>ode23s<\/tt> in evaluating the partial derivatives in each step. <tt>ode23tb<\/tt> also solves the problem with the fewest number of steps, outperforming <tt>ode15s<\/tt>. This problem is a good example of a stiff problem with a crude tolerance where <tt>ode23s<\/tt> and <tt>ode23tb<\/tt> can out perform <tt>ode15s<\/tt>. But practically speaking, all of the stiff solvers perform well on this problem and offer significant time savings when compared to <tt>ode45<\/tt>.<\/p><h4>Summary<a name=\"15751569-d461-4f8d-afa6-7cd24363a220\"><\/a><\/h4><p>Although all of the ODE solvers are capable of working on the same problems, it's recommended that you start with <tt>ode45<\/tt>. Then, if the problem exhibits signs of stiffness, <tt>ode15s<\/tt> is a good second choice. The other solvers then offer further refinement based on the properties of the specific problem and whether extra information (such as the Jacobian) can be provided.<\/p><h4>Comments<a name=\"afce60fc-d932-421c-b7f2-48491eba1b82\"><\/a><\/h4><p>Does your work involve the use of MATLAB's ODE solvers? If so, share your experience <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=1220#respond\">here<\/a>.<\/p><h4>Further Reading<a name=\"e197fc18-9fc4-4893-8c75-1316a54d1467\"><\/a><\/h4><p>[1] C. Moler, <a title=\"https:\/\/www.mathworks.com\/moler\/odes.pdf (link no longer works)\">Ordinary Differential Equations<\/a> <i>Numerical Computing with MATLAB<\/i> Electronic edition: The MathWorks, Inc., Natick, MA, 2004<\/p><p>[2] Shampine, L. F. and M. W. Reichelt, <a href=\"https:\/\/www.mathworks.com\/help\/pdf_doc\/otherdocs\/ode_suite.pdf\">The MATLAB ODE Suite<\/a> SIAM Journal on Scientific Computing, Vol. 18, 1997.<\/p><p>[3] C. Moler, <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2014\/06\/09\/ordinary-differential-equations-stiffness\/\">Ordinary Differential Equations: Stiffness<\/a> Cleve&#8217;s Corner: Cleve Moler on Mathematics and Computing, 2014<\/p><p>[4] L. F. Shampine, <i>Numerical Solution of Ordinary Differential Equations<\/i>, Chapman &amp; Hall, 1994<\/p><p>[5] Shampine, L. F., Gladwell, I. and S. Thompson, <i>Solving ODEs with MATLAB<\/i>, Cambridge University Press, 2003<\/p><p>[6] J. D. Lambert, <i>Numerical Methods for Ordinary Differential Systems<\/i>, New York: Wiley, 1992<\/p><p><i>Copyright 2015 The MathWorks, Inc.<\/i><\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_571dae2dd5d04d4fa56973dea75289df() {\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='571dae2dd5d04d4fa56973dea75289df ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 571dae2dd5d04d4fa56973dea75289df';\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 2015 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<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>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_571dae2dd5d04d4fa56973dea75289df()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2015a<br><\/p><\/div><!--\r\n571dae2dd5d04d4fa56973dea75289df ##### SOURCE BEGIN #####\r\n%% ODE Solver Selection in MATLAB\r\n% Today, I'd like to welcome Josh Meyer as this week's guest blogger. Josh\r\n% works on the Documentation team here at MathWorks, where he writes and\r\n% maintains some of the MATLAB Mathematics documentation. In this post,\r\n% Josh provides a bit of advice on how to choose which ODE solver to use.\r\n% Over to you, Josh...\r\n\r\n%% Initial Value Problems\r\n% There are 7 ordinary differential equation initial value problem solvers\r\n% in MATLAB:\r\n% \r\n% * |ode45|\r\n% * |ode23|\r\n% * |ode113|\r\n% * |ode15s|\r\n% * |ode23s|\r\n% * |ode23t|\r\n% * |ode23tb|\r\n% \r\n% _(note that |ode15i| is left out of this discussion because it solves its\r\n% own class of initial value problems: fully implicit ODEs of the form\r\n% $f(t,y,y') = 0$)_\r\n%\r\n% To choose between the solvers, it's first necessary to understand why one\r\n% solver might be better than another for a given problem.\r\n% \r\n% The ODE solvers in MATLAB all work on initial value problems of the form,\r\n% \r\n% $$y' = f \\left( t,y \\right)$$\r\n% \r\n% where $y' = dy\/dt$. There is also a more general form,\r\n% \r\n% $$ M(t,y) y' = f \\left( t,y \\right)$$\r\n% \r\n% where $M(t,y)$ is referred to as the _mass matrix_.\r\n%\r\n% Starting with the initial conditions $y_0$, and a period of time over\r\n% which the answer is to be obtained $(t_0,t_f)$, the solution is obtained\r\n% iteratively by using the results of previous steps according to the\r\n% solver's algorithm. At the first such step, the initial conditions\r\n% provide the necessary information that allows the integration to proceed.\r\n% The final result is that the ODE solver returns a vector of time steps\r\n% $t_0,t_1,...,t_f$ as well as the corresponding solution at each time step\r\n% $y_0,y_1,...,y_f$.\r\n% \r\n% Theoretically, this numerical solution technique is possible because of\r\n% the connection between differential equations and integrals provided by\r\n% the fundamental theorem of calculus:\r\n% \r\n% $$y(t + h) = y(t) + \\int_t^{t+h} f\\left( s,y(s) \\right) ds$$\r\n% \r\n% The problem of calculating $y(t+h)$ becomes a question of how to\r\n% approximate the integral on the right hand side. This is where different\r\n% solvers come in. Each different solver evaluates the integral using\r\n% different numerical techniques, and each solver makes trade-offs between\r\n% efficiency and accuracy.\r\n\r\n%% Example: Euler's Method\r\n% Euler's method is a simple ODE solver, but it provides an illustration of\r\n% the trade-offs between efficiency and accuracy in an ODE solver\r\n% algorithm. Suppose you want to solve\r\n% \r\n% $$y' = f(t,y) = 2t$$\r\n% \r\n% over the time span $[0,3]$ using the initial condition $y_0 = 0$. Each \r\n% step of Euler's method is computed with\r\n% \r\n% $$\\begin{array}{cl} y_{n+1} &= y_n + h f \\left(t_n,y_n\\right)\\\\ t_{n+1}&= \r\n% t_n + h\\end{array}$$\r\n% \r\n% Using $h=1$, the solution requires just three steps:\r\n% \r\n% $$\\begin{array}{cl} y_1 &= y_0 + f\\left(t_0,y_0\\right)=0\\\\y_2 &= y_1 +\r\n% f\\left(t_1,y_1 \\right)=2\\\\y_3 &= y_2 +  f\\left(t_2,y_2 \\right)=6\r\n% \\end{array}$$\r\n% \r\n% ... But is it accurate? \r\n% \r\n% Not really. The exact solution to this equation is \r\n% \r\n% $$y(t) = t^2$$\r\n% \r\n% Reducing the step size $h$ can improve the accuracy of the answer a bit,\r\n% but it also requires more steps to achieve the solution. To see this, the\r\n% below code solves this problem using Euler's method and compares the\r\n% answer to the analytic solution for several different values of $h$.\r\nclear, clc\r\nh = 1;\r\ntspan = [0 3];\r\nf = @(t,y) 2*t;\r\ndydt(1) = 0;\r\nt(1) = 0;\r\n\r\ny = @(t) t.^2;\r\nx = linspace(0,tspan(end));\r\nfigure\r\nplot(x,y(x))\r\nxlabel('t'), ylabel('y(t)')\r\nhold on\r\n\r\nwhile h > 0.1\r\n    for k = 2:tspan(end)\/h+1\r\n        dydt(k) = dydt(k-1) + h*f(t(k-1),dydt(k-1));\r\n        t(k) = t(k-1) + h;\r\n    end\r\n    \r\n    plot(t,dydt,'-o')\r\n    h = h\/2;\r\nend\r\nlegend('Exact Solution', 'h=1','h=0.5','h=0.25','h=0.125',...\r\n    'Location','NorthWest')\r\ntitle('Solution of y''=2t using Euler''s method with several step sizes')\r\n\r\n\r\n%% Improving on Euler's Method\r\n% Using smaller and smaller step sizes turns out to not be a good idea,\r\n% since the algorithm loses efficiency. For any reasonable problem such a\r\n% solver would be very slow. Also, Euler's method has a few inherent\r\n% problems. Since the slope of $y$ is evaluated only once at the beginning\r\n% of each interval, this solver only produces exact answers for constant\r\n% functions. There is also no way to estimate the error, so the solver\r\n% needs to use fixed step sizes.\r\n%\r\n% So, one way to improve on Euler's method is to evaluate $y'$ more often\r\n% in each step. This provides intermediate slopes that give a better idea\r\n% of what the function is doing within each interval, allowing the solver\r\n% to produce exact answers for higher order problems. For example, if you\r\n% add an evaluation of the slope halfway across each interval to Euler's\r\n% method, then the result is called the _midpoint rule_, which produces\r\n% exact integrations for linear functions:\r\n%\r\n% $$\\begin{array}{cl} s_1 &= f(t_n,y_n)\\\\ s_2 &= f\\left( t_n + \\frac{h}{2},\r\n% y_n + \\frac{h}{2}s_1\\right)\\\\y_{n+1} &= y_n + hs_2\\\\t_{n+1} &=\r\n% t_n+h\\end{array}$$\r\n%\r\n% If you evaluate the slope four times in each interval, you get the\r\n% _classical Runge-Kutta_ algorithm (a.k.a. _RK4_), which is a piece of the\r\n% |ode45| algorithm. This algorithm produces exact integrations for cubic\r\n% functions (and if $f$ is only a function of $t$, then $s_2=s_3$ and this\r\n% is the same as _Simpson's rule_ for quadrature):\r\n%\r\n% $$\\begin{array}{cl} s_1 &= f(t_n,y_n)\\\\s_2 &=\r\n% f\\left(t_n+\\frac{h}{2},y_n+\\frac{h}{2}s_1\\right)\\\\s_3 &=\r\n% f\\left(t_n+\\frac{h}{2},y_n+\\frac{h}{2}s_2\\right)\\\\ s_4 &=\r\n% f\\left(t_n+h,y_n+hs_3\\right)\\\\y_{n+1} &= y_n +\r\n% \\frac{h}{6}\\left(s_1+2s_2+2s_3+s_4\\right)\\\\t_{n+1} &= t_n +\r\n% h\\end{array}$$\r\n%\r\n% Runge-Kutta algorithms are all _single-step_ solvers, since each step\r\n% only depends on the result of the previous step. |ode45|, |ode23|,\r\n% |ode23s|, |ode23t|, and |ode23tb| all employ single-step algorithms.\r\n% Multi-step algorithms, such as those employed by |ode113| and |ode15s|,\r\n% use the results of several past steps.\r\n% \r\n% Sophisticated ODE solvers, like the ones in MATLAB, also estimate the\r\n% error in each step to determine how big the next step size should be.\r\n% This is another improvement over the fixed step sizes used above, since a\r\n% solver that does more work per step is able to compensate by taking steps\r\n% of varying size. The error estimate used to determine the step size is\r\n% typically obtained by comparing the results of two different methods.\r\n% MATLAB's ODE solvers follow a naming convention that reveals information\r\n% about which methods they use. |ode45| compares the results of a 4th-order\r\n% Runge-Kutta method and a 5th-order Runge-Kutta method to determine the\r\n% error. Similarly, |ode23| uses a 2nd-order and 3rd-order Runge-Kutta\r\n% comparison. So, in general, the smaller the number |odeNN|, the looser\r\n% the solver's error tolerance is.\r\n%\r\n% It should be no surprise, then, that |ode45| obtains a very accurate\r\n% answer for the equation we solved before with Euler's method. |ode45|\r\n% is MATLAB's general purpose ODE solver, and it is the first solver you\r\n% should use for most problems.\r\ny = @(t) t.^2;\r\nx = linspace(0,3);\r\nfigure\r\nplot(x,y(x))\r\nxlabel('t'), ylabel('y(t)')\r\nhold on\r\n[t,y] = ode45(@(t,y) 2*t, [0 3], 0);\r\nplot(t,y,'o')\r\nxlabel('t'), ylabel('y(t)')\r\ntitle('Solution of y''=2t using ode45')\r\n\r\n%% Stiff Differential Equations\r\n% For some ODE problems, the step size taken by the solver is forced down\r\n% to an unreasonably small level in comparison to the interval of\r\n% integration, even in a region where the solution curve is smooth. These\r\n% step sizes can be so small that traversing a short time interval might\r\n% require millions of evaluations. This can lead to the solver failing the\r\n% integration, but even if it succeeds it will take a very long time to do\r\n% so. \r\n%\r\n% Equations that cause this behavior in ODE solvers are said to be _stiff_.\r\n% This is a nod to the fact that the equations are stubborn and not easily\r\n% evaluated with numerical techniques. The problem that stiff ODEs pose is\r\n% that explicit solvers (such as |ode45|) are untenably slow in achieving a\r\n% solution. This is why |ode45| is classified as a nonstiff solver along\r\n% with |ode23| and |ode113|. These solvers all struggle to integrate stiff\r\n% equations.\r\n%\r\n% Equation stiffness resists a precise definition, because there are\r\n% several factors that cause it. Stiffness results from a combination of\r\n% the specific equations, the ODE solver being used, the initial\r\n% conditions, and the error tolerance used by the solver. The following\r\n% statements about stiffness, attributed to Lambert [6], are exhibited by\r\n% many examples of stiff ODEs, but counterexamples also exist, so they are\r\n% not true definitions of stiffness:\r\n% \r\n% # A linear constant coefficient system is stiff if all of its\r\n% eigenvalues have negative real part and the stiffness ratio [of the\r\n% largest and smallest eigenvalues] is large.\r\n% # Stiffness occurs when the mathematical problem is stable, and yet\r\n% stability requirements, rather than those of accuracy, severely constrain\r\n% the step length.\r\n% # Stiffness occurs when some components of the solution decay much more\r\n% rapidly than others.\r\n%\r\n% A common theme among these statements is that stiffness can result from a\r\n% difference in scaling somewhere in the problem. This difference in scale\r\n% (for example, if the Jacobian $J = \\partial f_n\/\\partial y_i$ has a large\r\n% ratio of negative eigenvalues) constrains the step size that the solver\r\n% can take in performing the integration. Tiny step sizes become necessary\r\n% in order to preserve any notion of error tolerance or stability in the\r\n% solution.\r\n% \r\n% For example, equations describing chemical reactions frequently display\r\n% stiffness, since it is common for components of the solution to vary on\r\n% drastically different time scales (reactions occurring at the same time\r\n% that are both very slow and very fast).\r\n% \r\n% However, there are solvers specifically designed to work on stiff ODEs.\r\n% Solvers that are designed for stiff problems typically do more work per\r\n% step, and the pay-off is that they are able to take much larger steps and\r\n% enjoy improved numerical stability compared to the nonstiff solvers.\r\n% Stiff solvers are implicit, because the computation of $y$ requires the\r\n% use of linear algebra to solve systems of linear equations. The Jacobian\r\n% is used to estimate the local behavior of the ODE as the integration\r\n% proceeds, so supplying the analytical Jacobian can improve the\r\n% performance of MATLAB's stiff ODE solvers.\r\n%\r\n% This is just a cursory treatment of stiffness, because it is a complex\r\n% topic. See\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2014\/06\/09\/ordinary-differential-equations-stiffness\/\r\n% Ordinary Differential Equations: Stiffness> for a more in-depth look.\r\n%\r\n% To summarize, the nonstiff solvers in MATLAB are:\r\n%\r\n% * |ode45|\r\n% * |ode23|\r\n% * |ode113|\r\n%\r\n% The stiff solvers are (when |ode45| is slow):\r\n%\r\n% * |ode15s|\r\n% * |ode23s|\r\n% * |ode23t|\r\n% * |ode23tb|\r\n%\r\n% It should be noted that nonstiff solvers do _work_ on stiff problems, it\r\n% is just that they are exceptionally slow. Similarly, solvers designed for\r\n% stiff problems can work on nonstiff problems, but since they do more work\r\n% per step they are less efficient than their nonstiff counterparts when\r\n% that extra work isn't necessary. So equation stiffness is a matter of\r\n% solver efficiency, and the goal is to strike the right balance between\r\n% accuracy of the solution and work done in each step by the solver.\r\n\r\n%% Solver Recommendations\r\n% The following recommendations are adapted from the\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/math\/choose-an-ode-solver.html\r\n% MATLAB Mathematics documentation> :\r\n%\r\n% * |ode45| is MATLAB's general purpose single-step ODE solver. This should\r\n%   be the first solver you use for most problems.\r\n% \r\n% For *nonstiff* problems:\r\n% \r\n% * |ode23| is another single-step solver that can be more efficient than\r\n%   |ode45| if the problem permits a crude error tolerance. This looser\r\n%   error tolerance can also accommodate some mildly stiff problems.\r\n% * |ode113| is a multi-step solver, and is preferred over |ode45| if the\r\n%   function is expensive to evaluate, or for smooth problems where high\r\n%   precision is required. For example, |ode113| excels with orbital\r\n%   dynamics and celestial mechanics problems.\r\n% \r\n% For *stiff* problems (where |ode45| is slow):\r\n% \r\n% * |ode15s| is a multi-step solver that is MATLAB's general purpose solver\r\n%   for stiff problems. Use |ode15s| if |ode45| fails or struggles to\r\n%   complete the integration in a reasonable amount of time. |ode15s| is\r\n%   also the primary solver for DAEs, which are identified as ODEs with a\r\n%   singular mass matrix.\r\n% * For stiff problems with crude error tolerances, |ode23s|, |ode23t|, and\r\n%   |ode23tb| provide more efficient alternatives to |ode15s| since they\r\n%   are single-step solvers. The efficiency of |ode23s| can be\r\n%   significantly improved by providing the Jacobian, since |ode23s|\r\n%   evaluates the Jacobian in each step.\r\n% * |ode23s| only works on ODEs with a mass matrix if the mass matrix is\r\n%   constant (not time- or state-dependent).\r\n% * |ode15s| and |ode23t| are the only solvers that solve DAEs of index 1.\r\n%   \r\n%\r\n% Here is a graphic that captures the basic recommendations. In most cases,\r\n% the only choice in solver you will need to make is to use |ode15s|\r\n% instead of |ode45|. \r\n%\r\n% <<selector.png>>\r\n%\r\n\r\n\r\n%% Example 1: Damped Pendulum\r\n% The equation of motion for a damped pendulum is,\r\n%\r\n% $$\\ddot{\\theta} = -\\frac{b}{m}\\dot{\\theta}-\\frac{mg}{L\\left(m-2b\\right)}\r\n% \\sin \\theta$$\r\n%\r\n% where $g$ is the gravitational constant, $m$ the mass of the bob, $L$ the\r\n% length of the string, and $b$ is a damping coefficient. The goal is to\r\n% solve for $\\theta$, the angle that the pendulum deviates from the\r\n% vertical, and $\\theta '$, the rate at which the angle changes.\r\n% \r\n% Some natural initial conditions would be $\\theta_0 = \\pi\/4$ and $\\theta\r\n% '_0 = 0$, indicating that you lift the pendulum up to a 45 degree angle\r\n% before letting go, and it has no initial angular velocity. Due to the\r\n% damping coefficient, you would expect the pendulum to slowly lose\r\n% momentum and go back down to rest.\r\n%\r\n% The file |pendulumODE.m| reformulates the problem as a coupled system of\r\n% first-order ODEs:\r\n%\r\n% $$\\begin{array}{cl} y_1' &= y_2\\\\y_2' &= -\\frac{b}{m} y_2\r\n% -\\frac{mg}{L(m-2b)}sin(y_1)\\end{array}$$\r\n% \r\n% then solves using |ode45|, |ode15s|, |ode23|, and |ode113|. The solutions\r\n% for $y_1 = \\theta$ are plotted, and the file returns the stats for each\r\n% solver. As is always the case when displaying execution times, \"the\r\n% timings displayed can vary\".\r\n%\r\n% <include>pendulumODE.m<\/include>\r\n%\r\npendulumODE\r\n\r\n%%\r\n% The solvers all perform well, but the damped pendulum is a good example\r\n% of a nonstiff problem where |ode45| performs nicely. In this case\r\n% |ode15s| needs to do extra work in order to achieve an inferior solution.\r\n\r\n%% Example 2: van der Pol Oscillator\r\n% The van der Pol Oscillator equation becomes stiff in certain intervals\r\n% when the nonlinear parameter $\\mu$ is large:\r\n%\r\n% $$\\ddot{y} - \\mu \\left(1-y^2\\right)\\dot{y}+y=0$$\r\n%\r\n% The nonlinearity of this equation is contained entirely in the term that\r\n% involves $\\mu$: notice that if $\\mu=0$, the equation reduces to that of a\r\n% simple harmonic oscillator, which has regular periodic behavior.\r\n%\r\n% Attempting to solve this equation using |ode45| is met with severe\r\n% resistance, requiring millions of evaluations and 30+ minutes of\r\n% execution (I stopped execution after 35 minutes). Since the problem is\r\n% clearly stiff, this example compares the stiff solvers.\r\n%\r\n% The file |vanderpolODE.m| finds the solution for $\\mu=1000$ using\r\n% |ode15s|, |ode23s|, |ode23t|, and |ode23tb|. The function file\r\n% |vdp1000.m| ships with MATLAB and encodes this equation as a coupled\r\n% system of first-order ODEs:\r\n%\r\n% $$\\begin{array}{cl}y_1' &= y_2\\\\y_2' &= \\mu \\left( 1-y_1^2 \\right)y_2 -\r\n% y_1\\end{array}$$\r\n%\r\n% The Jacobian is supplied to assist the solvers, and its use is reflected\r\n% in the number of partial derivative evaluations.\r\n%\r\n% <include>vanderpolODE.m<\/include>\r\n%\r\nvanderpolODE\r\n\r\n%%\r\n% The plots are of the solutions for $y_1$. For this problem, |ode23s|\r\n% executes quickest and with the least number of failed steps. The supplied\r\n% Jacobian greatly assists |ode23s| in evaluating the partial derivatives\r\n% in each step. |ode23tb| also solves the problem with the fewest number of\r\n% steps, outperforming |ode15s|. This problem is a good example of a stiff\r\n% problem with a crude tolerance where |ode23s| and |ode23tb| can out\r\n% perform |ode15s|. But practically speaking, all of the stiff solvers\r\n% perform well on this problem and offer significant time savings when\r\n% compared to |ode45|.\r\n\r\n%% Summary\r\n% Although all of the ODE solvers are capable of working on the same\r\n% problems, it's recommended that you start with |ode45|. Then, if the\r\n% problem exhibits signs of stiffness, |ode15s| is a good second choice.\r\n% The other solvers then offer further refinement based on the properties\r\n% of the specific problem and whether extra information (such as the\r\n% Jacobian) can be provided.\r\n\r\n%% Comments\r\n% Does your work involve the use of MATLAB's ODE solvers? If so, share your\r\n% experience <https:\/\/blogs.mathworks.com\/loren\/?p=1220#respond here>.\r\n\r\n%% Further Reading\r\n% [1] C. Moler, <https:\/\/www.mathworks.com\/moler\/odes.pdf Ordinary Differential\r\n% Equations> _Numerical Computing with MATLAB_ Electronic edition: The\r\n% MathWorks, Inc., Natick, MA, 2004\r\n%\r\n% [2] Shampine, L. F. and M. W. Reichelt,\r\n% <https:\/\/www.mathworks.com\/help\/pdf_doc\/otherdocs\/ode_suite.pdf The\r\n% MATLAB ODE Suite> SIAM Journal on Scientific Computing, Vol. 18, 1997.\r\n%\r\n% [3] C. Moler,\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2014\/06\/09\/ordinary-differential-equations-stiffness\/\r\n% Ordinary Differential Equations: Stiffness> Cleve\u00e2\u20ac\u2122s Corner: Cleve Moler\r\n% on Mathematics and Computing, 2014\r\n%\r\n% [4] L. F. Shampine, _Numerical Solution of Ordinary Differential\r\n% Equations_, Chapman & Hall, 1994\r\n%\r\n% [5] Shampine, L. F., Gladwell, I. and S. Thompson, _Solving ODEs with\r\n% MATLAB_, Cambridge University Press, 2003\r\n%\r\n% [6] J. D. Lambert, _Numerical Methods for Ordinary Differential Systems_,\r\n% New York: Wiley, 1992\r\n\r\n%%\r\n% _Copyright 2015 The MathWorks, Inc._\r\n##### SOURCE END ##### 571dae2dd5d04d4fa56973dea75289df\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/OrdinaryDifferentialEquationsInMATLAB_04.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>Today, I'd like to welcome Josh Meyer as this week's guest blogger. Josh works on the Documentation team here at MathWorks, where he writes and maintains some of the MATLAB Mathematics documentation. In this post, Josh provides a bit of advice on how to choose which ODE solver to use. Over to you, Josh...... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2015\/09\/23\/ode-solver-selection-in-matlab\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[26,31,1],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1220"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=1220"}],"version-history":[{"count":3,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1220\/revisions"}],"predecessor-version":[{"id":1223,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1220\/revisions\/1223"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=1220"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=1220"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=1220"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}