09427372e89e81f6f5350fa3d09dc093a3c7f198
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) Scilab Enterprises - 2012 - Paul Bignier
5  *
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.
9  * The terms are also available at
10  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
11  -->
13     <refnamediv>
14         <refname>Implicit Runge-Kutta 4(5)</refname>
15         <refpurpose>
16             <emphasis>Implicit Runge-Kutta</emphasis> is a numerical solver providing an efficient and stable implicit method to solve Ordinary Differential Equations (ODEs) Initial Value Problems. Called by <link linkend="xcos">xcos</link>.
17         </refpurpose>
18     </refnamediv>
19     <refsection>
20         <title>Description</title>
21         <para>
22             <emphasis>Runge-Kutta</emphasis> is a numerical solver providing an efficient and stable fixed-size step method to solve Initial Value Problems of the form:
23         </para>
24         <para>
25             <latex>
26                 \begin{eqnarray}
27                 \dot{y} = f(t,y), \hspace{3 mm} y(t_0) = y_0, \hspace{3 mm} y \in R^N
28                 \end{eqnarray}
29             </latex>
30         </para>
31         <para>
32             <emphasis>CVode</emphasis> and <emphasis>IDA</emphasis> use variable-size steps for the integration.
33         </para>
34         <para>
35             A drawback of that is the unpredictable computation time. With <emphasis>Runge-Kutta</emphasis>, we do not adapt to the complexity of the problem, but we guarantee a stable computation time.
36         </para>
37         <para>
38             This method being implicit, it can be used on stiff problems.
39         </para>
40         <para>
41             It is an enhancement of the backward Euler method, which approximates
42             <emphasis>
43                 y<subscript>n+1</subscript>
44             </emphasis>
45             by computing
46             <emphasis>
47                 f(t<subscript>n</subscript>+h, y<subscript>n+1</subscript>)
48             </emphasis>
49             and truncating the Taylor expansion.
50         </para>
51         <para>
52             The implemented scheme is inspired from the "Low-Dispersion Low-Dissipation Implicit Runge-Kutta Scheme" (see bottom for link).
53         </para>
54         <para>
56         </para>
57         <para>
58             An important difference of <emphasis>implicit Runge-Kutta</emphasis> with the previous methods is that it computes up to the fourth derivative of <emphasis>y</emphasis>, while the others mainly use linear combinations of <emphasis>y</emphasis> and <emphasis>y'</emphasis>.
59         </para>
60         <para>
61             Here, the next value is determined by the present value
62             <emphasis>
63                 y<subscript>n</subscript>
64             </emphasis>
65             plus the weighted average of three increments, where each increment is the product of the size of the interval, <emphasis>h</emphasis>, and an estimated slope specified by the function <emphasis>f(t,y)</emphasis>. They are distributed approximately equally on the interval.
66             <itemizedlist>
67                 <listitem>
68                     <emphasis>k1</emphasis> is the increment based on the slope near the quarter of the interval, using
69                     <emphasis>
70                         y<subscript>n</subscript>+ a11*h*k1,
71                     </emphasis>
72                     ,
73                 </listitem>
74                 <listitem>
75                     <emphasis>k2</emphasis> is the increment based on the slope near the midpoint of the interval, using
76                     <emphasis>
77                         y<subscript>n</subscript> + a21*h*k1 + a22*h*k2,
78                     </emphasis>
79                     ,
80                 </listitem>
81                 <listitem>
82                     <emphasis>k3</emphasis> is the increment based on the slope near the third quarter of the interval, using
83                     <emphasis>
84                         y<subscript>n</subscript> + a31*h*k1 + a32*h*k2 + a33*h*k3.
85                     </emphasis>
86                 </listitem>
87             </itemizedlist>
88         </para>
89         <para>
90             We see that the computation of <emphasis>ki</emphasis> requires <emphasis>ki</emphasis>, thus necessitating the use of a nonlinear solver (here, fixed-point iterations).
91         </para>
92         <para>
93             First, we set
94             <emphasis>
95                 k0 = h * f(t<subscript>n</subscript>, y<subscript>n</subscript>)
96             </emphasis>
97             as first guess for all the <emphasis>ki</emphasis>, to get updated <emphasis>ki</emphasis> and a first value for
98             <emphasis>
99                 y<subscript>n+1</subscript>
100             </emphasis>
101             .
102         </para>
103         <para>
104             Next, we save and recompute
105             <emphasis>
106                 y<subscript>n+1</subscript>
107             </emphasis>
108             with those new <emphasis>ki</emphasis>.
109         </para>
110         <para>
111             Then, we compare the two
112             <emphasis>
113                 y<subscript>n+1</subscript>
114             </emphasis>
115             and recompute it until its difference with the last computed one is inferior to the simulation parameter <emphasis>reltol</emphasis>.
116         </para>
117         <para>
118             This process adds a significant computation time to the method, but greatly improves stability.
119         </para>
120         <para>
121             We can see that with the <emphasis>ki</emphasis>, we progress in the derivatives of
122             <emphasis>
123                 y<subscript>n</subscript>
124             </emphasis>
125             . So in <emphasis>k3</emphasis>, we are approximating
126             <emphasis>
127                 y<superscript>(3)</superscript><subscript>n</subscript>
128             </emphasis>
129             , thus making an error in
130             <emphasis>
131                 O(h<superscript>4</superscript>)
132             </emphasis>
133             . But choosing the right coefficients in the computation of the <emphasis>ki</emphasis> (notably the
134             <emphasis>
135                 a<subscript>ij</subscript>
136             </emphasis>
137             ) makes us gain an order, thus making a per-step total error in
138             <emphasis>
139                 O(h<superscript>5</superscript>)
140             </emphasis>
141             .
142         </para>
143         <para>
144             So the total error is
145             <emphasis>
146                 number of steps * O(h<superscript>5</superscript>)
147             </emphasis>
148             . And since <emphasis>number of steps = interval size / h</emphasis> by definition, the total error is in
149             <emphasis>
150                 O(h<superscript>4</superscript>)
151             </emphasis>
152             .
153         </para>
154         <para>
155             That error analysis baptized the method <emphasis>implicit Runge-Kutta 4(5)</emphasis>:
156             <emphasis>
157                 O(h<superscript>5</superscript>)
158             </emphasis>
159             per step,
160             <emphasis>
161                 O(h<superscript>4</superscript>)
162             </emphasis>
163             in total.
164         </para>
165         <para>
167             <emphasis>
168                 10<superscript>-3</superscript>
169             </emphasis>
170             , rounding errors sometimes come into play as we approach
171             <emphasis>
172                 4*10<superscript>-4</superscript>
173             </emphasis>
174             . Indeed, the interval splitting cannot be done properly and we get capricious results.
175         </para>
176     </refsection>
177     <refsection>
178         <title>Examples</title>
179         <para>
181                 <inlinemediaobject>
182                     <imageobject>
183                         <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.zcos" valign="middle"/>
184                     </imageobject>
185                 </inlinemediaobject>
187             <scilab:image><![CDATA[
190 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.zcos");
191 scs_m.props.tol(2) = 10^-5;
192 scs_m.props.tol(6) = 7;
193 scs_m.props.tol(7) = 10^-2;
194 try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
195 ]]></scilab:image>
196         </para>
197         <para>
198             The integral block returns its continuous state, we can evaluate it with <emphasis>implicit RK</emphasis> by running the example:
199         </para>
200         <para>
201             <programlisting language="example"><![CDATA[
202       // Import the diagram and set the ending time
205       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.zcos");
206       scs_m.props.tf = 5000;
208       // Select the solver implicit RK and set the precision
209       scs_m.props.tol(2) = 10^-10;
210       scs_m.props.tol(6) = 7;
211       scs_m.props.tol(7) = 10^-2;
213       // Start the timer, launch the simulation and display time
214       tic();
215       try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
216       t = toc();
217       disp(t, "Time for implicit Runge-Kutta:");
218       ]]></programlisting>
219         </para>
220         <para>
221             The Scilab console displays:
222             <screen><![CDATA[
223 Time for implicit Runge-Kutta:
224  8.911
225             ]]></screen>
226         </para>
227         <para>
228             Now, in the following script, we compare the time difference between <emphasis>implicit RK</emphasis> and <emphasis>CVode</emphasis> by running the example with the five solvers in turn:
230                 Open the script
232         </para>
233         <para>
234             <screen><![CDATA[
235 Time for BDF / Newton:
236  18.894
238 Time for BDF / Functional:
239  18.382
241 Time for Adams / Newton:
242  10.368
244 Time for Adams / Functional:
245  9.815
247 Time for implicit Runge-Kutta:
248  6.652
249             ]]></screen>
250         </para>
251         <para>
252             These results show that on a nonstiff problem, for relatively same precision required and forcing the same step size, <emphasis>implicit Runge-Kutta</emphasis> competes with <emphasis>Adams / Functional</emphasis>.
253         </para>
254         <para>
255             Variable-size step ODE solvers are not appropriate for deterministic real-time applications because the computational overhead of taking a time step varies over the course of an application.
256         </para>
257     </refsection>
258     <refsection>
259         <title>See Also</title>
260         <simplelist type="inline">
261             <member>
263             </member>
264             <member>
266             </member>
267             <member>
269             </member>
270             <member>
272             </member>
273             <member>
275             </member>
276             <member>
278             </member>
279             <member>
281             </member>
282             <member>
284             </member>
285             <member>
287             </member>
288             <member>
290             </member>
291             <member>
293             </member>
294         </simplelist>
295     </refsection>
296     <refsection>
297         <title>Bibliography</title>
298         <para>
299             Journal of Computational Physics, Volume 233, January 2013, Pages 315-323 <ulink url="http://dl.acm.org/citation.cfm?id=2397727.2397975&amp;coll=DL&amp;dl=GUIDE&amp;CFID=295003167&amp;CFTOKEN=37159263">A low-dispersion and low-dissipation implicit Runge-Kutta scheme</ulink>
300         </para>
301         <para>