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-en.txt
11  -->
13     <refnamediv>
14         <refname>Runge-Kutta 4(5)</refname>
15         <refpurpose>
16             <emphasis>Runge-Kutta</emphasis> is a numerical solver providing an efficient explicit method to solve Ordinary Differential Equations (ODEs) Initial Value Problems.
17         </refpurpose>
18     </refnamediv>
19     <refsection>
20         <title>Description</title>
21         <para>
22             Called by <link linkend="xcos">xcos</link>, <emphasis>Runge-Kutta</emphasis> is a numerical solver providing an efficient 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 Runge-Kutta, we do not adapt to the complexity of the problem, but we guarantee a stable computation time.
36         </para>
37         <para>
38             As of now, this method is explicit, so it is not concerned with Newton or Functional iterations, and not advised for stiff problems.
39         </para>
40         <para>
41             It is an enhancement of the Euler method, which approximates
42             <emphasis>
43                 y<subscript>n+1</subscript>
44             </emphasis>
45             by truncating the Taylor expansion.
46         </para>
47         <para>
49         </para>
50         <para>
51             An important difference of Runge-Kutta with the previous methods is that it computes up to the fourth derivative of <emphasis>y</emphasis>, while the others only use linear combinations of <emphasis>y</emphasis> and <emphasis>y'</emphasis>.
52         </para>
53         <para>
54             Here, the next value is determined by the present value
55             <emphasis>
56                 y<subscript>n</subscript>
57             </emphasis>
58             plus the weighted average of four 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> :
59             <itemizedlist>
60                 <listitem>
61                     <emphasis>k1</emphasis> is the increment based on the slope at the beginning of the interval, using
62                     <emphasis>
63                         y<subscript>n</subscript>
64                     </emphasis>
65                     (Euler's method),
66                 </listitem>
67                 <listitem>
68                     <emphasis>k2</emphasis> is the increment based on the slope at the midpoint of the interval, using
69                     <emphasis>
70                         y<subscript>n</subscript> + h*k1/2
71                     </emphasis>
72                     ,
73                 </listitem>
74                 <listitem>
75                     <emphasis>k3</emphasis> is again the increment based on the slope at the midpoint, but now using
76                     <emphasis>
77                         y<subscript>n</subscript> + h*k2/2
78                     </emphasis>
79                 </listitem>
80                 <listitem>
81                     <emphasis>k4</emphasis> is the increment based on the slope at the end of the interval, using
82                     <emphasis>
83                         y<subscript>n</subscript> + h*k3
84                     </emphasis>
85                 </listitem>
86             </itemizedlist>
87         </para>
88         <para>
89             We can see that with the <emphasis>ki</emphasis>, we progress in the derivatives of
90             <emphasis>
91                 y<subscript>n</subscript>
92             </emphasis>
93             . So in <emphasis>k4</emphasis>, we are approximating
94             <emphasis>
95                 y<superscript>(4)</superscript><subscript>n</subscript>
96             </emphasis>
97             , thus making an error in
98             <emphasis>
99                 O(h<superscript>5</superscript>)
100             </emphasis>
101             .
102         </para>
103         <para>
104             So the total error is
105             <emphasis>
106                 number of steps * O(h<superscript>5</superscript>)
107             </emphasis>
108             . And since <emphasis>number of steps = interval size / h</emphasis> by definition, the total error is in
109             <emphasis>
110                 O(h<superscript>4</superscript>)
111             </emphasis>
112             .
113         </para>
114         <para>
115             That error analysis baptized the method <emphasis>Runge-Kutta 4(5)</emphasis>,
116             <emphasis>
117                 O(h<superscript>5</superscript>)
118             </emphasis>
119             per step,
120             <emphasis>
121                 O(h<superscript>4</superscript>)
122             </emphasis>
123             in total.
124         </para>
125         <para>
127             <emphasis>
128                 10<superscript>-3</superscript>
129             </emphasis>
130             , rounding errors sometimes come into play as we approach <emphasis>
131                 4*10<superscript>-4</superscript>
132             </emphasis>
133             . Indeed, the interval splitting cannot be done properly and we get capricious results.
134         </para>
135     </refsection>
136     <refsection>
137         <title>Examples</title>
138         <para>
140                 <inlinemediaobject>
141                     <imageobject>
142                         <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.xcos" valign="middle"/>
143                     </imageobject>
144                 </inlinemediaobject>
146             <scilab:image><![CDATA[
149 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
150 scs_m.props.tol(6) = 6;
151 scs_m.props.tol(7) = 10^-2;
152 try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
153 ]]></scilab:image>
154         </para>
155         <para>
156             The integral block returns its continuous state, we can evaluate it with Runge-Kutta by running the example :
157         </para>
158         <para>
159             <programlisting language="example"><![CDATA[
160       // Import the diagram and set the ending time
163       importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
164       scs_m.props.tf = 5000;
166       // Select the solver Runge-Kutta and set the precision
167       scs_m.props.tol(6) = 6;
168       scs_m.props.tol(7) = 10^-2;
170       // Start the timer, launch the simulation and display time
171       tic();
172       try xcos_simulate(scs_m, 4); catch disp(lasterror()); end
173       t = toc();
174       disp(t, "Time for Runge-Kutta :");
175       ]]></programlisting>
176         </para>
177         <para>
178             The Scilab console displays :
179             <screen><![CDATA[
180 Time for Runge-Kutta :
181  15.624
182             ]]></screen>
183         </para>
184         <para>
185             Now, in the following script, we compare the time difference between Runge-Kutta and Sundials by running the example with the five solvers in turn :
187                 Open the script
189         </para>
190         <para>
191             <screen><![CDATA[
192 Time for BDF / Newton :
193  19.894
195 Time for BDF / Functional :
196  19.382
198 Time for Adams / Newton :
199  11.255
201 Time for Adams / Functional :
202  10.233
204 Time for Runge-Kutta :
205  5.502
206             ]]></screen>
207         </para>
208         <para>
209             These results show that on a nonstiff problem, for relatively same precision required and forcing the same step size, Runge-Kutta is faster.
210         </para>
211         <para>
212             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.
213         </para>
214     </refsection>
215     <refsection>
216         <title>See Also</title>
217         <simplelist type="inline">
218             <member>
220             </member>
221             <member>
223             </member>
224             <member>
226             </member>
227             <member>
229             </member>
230             <member>
232             </member>
233             <member>
235             </member>
236             <member>
238             </member>
239             <member>
241             </member>
242             <member>
244             </member>
245             <member>
247             </member>
248             <member>
250             </member>
251         </simplelist>
252     </refsection>
253     <refsection>
254         <title>Bibliography</title>
255         <para>