<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>Crystal Space 1.2.1: csgeom/odesolver.h Source File (Crystal Space 1.2.1 Public API Reference)</title>
<link href="tabs.css" rel="stylesheet" type="text/css">
<link href="doxygen.css" rel="stylesheet" type="text/css">
</head><body>
<table border="0" cellpadding="0" cellspacing="0" width="100%" class="head">
<tr height="59">
<td class="head" width="202" valign="bottom" style="padding-left:0;"><a href="http://www.crystalspace3d.org/"><img src="csblur.png" width="236" height="59" alt="CrystalSpace" border="0"></a></td>
<td class="head"><h2>Public API Reference</h2></td>
</tr>
<tr height="11">
<td colspan="2" class="headshadow" valign="top" style="padding-left:0;"><img src="csblurb.png" width="236" height="11" alt="" border="0"></td>
</tr>
</table>
<div class="content">
<!-- Generated by Doxygen 1.5.3 -->
<div class="tabs">
<ul>
<li><a href="index.html"><span>Main Page</span></a></li>
<li><a href="modules.html"><span>Modules</span></a></li>
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
<li><a href="classes.html"><span>Classes</span></a></li>
<li class="current"><a href="files.html"><span>Files</span></a></li>
<li><a href="pages.html"><span>Related Pages</span></a></li>
</ul>
</div>
<h1>csgeom/odesolver.h</h1><a href="odesolver_8h.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">/*</span>
<a name="l00002"></a>00002 <span class="comment"> Copyright (C) 2006 by Marten Svanfeldt</span>
<a name="l00003"></a>00003 <span class="comment"></span>
<a name="l00004"></a>00004 <span class="comment"> This library is free software; you can redistribute it and/or</span>
<a name="l00005"></a>00005 <span class="comment"> modify it under the terms of the GNU Library General Public</span>
<a name="l00006"></a>00006 <span class="comment"> License as published by the Free Software Foundation; either</span>
<a name="l00007"></a>00007 <span class="comment"> version 2 of the License, or (at your option) any later version.</span>
<a name="l00008"></a>00008 <span class="comment"></span>
<a name="l00009"></a>00009 <span class="comment"> This library is distributed in the hope that it will be useful,</span>
<a name="l00010"></a>00010 <span class="comment"> but WITHOUT ANY WARRANTY; without even the implied warranty of</span>
<a name="l00011"></a>00011 <span class="comment"> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU</span>
<a name="l00012"></a>00012 <span class="comment"> Library General Public License for more details.</span>
<a name="l00013"></a>00013 <span class="comment"></span>
<a name="l00014"></a>00014 <span class="comment"> You should have received a copy of the GNU Library General Public</span>
<a name="l00015"></a>00015 <span class="comment"> License along with this library; if not, write to the Free</span>
<a name="l00016"></a>00016 <span class="comment"> Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.</span>
<a name="l00017"></a>00017 <span class="comment">*/</span>
<a name="l00018"></a>00018
<a name="l00023"></a>00023 <span class="preprocessor">#ifndef __CS_CSGEOM_ODESOLVER_H__</span>
<a name="l00024"></a>00024 <span class="preprocessor"></span><span class="preprocessor">#define __CS_CSGEOM_ODESOLVER_H__</span>
<a name="l00025"></a>00025 <span class="preprocessor"></span>
<a name="l00026"></a>00026 <span class="preprocessor">#include "<a class="code" href="vector3_8h.html" title="3D vector.">csgeom/vector3.h</a>"</span>
<a name="l00027"></a>00027
<a name="l00028"></a><a class="code" href="namespaceCS.html">00028</a> <span class="keyword">namespace </span>CS
<a name="l00029"></a>00029 {
<a name="l00030"></a><a class="code" href="namespaceCS_1_1Math.html">00030</a> <span class="keyword">namespace </span>Math
<a name="l00031"></a>00031 {
<a name="l00032"></a>00032
<a name="l00046"></a><a class="code" href="classCS_1_1Math_1_1Ode45.html">00046</a> <span class="keyword">class </span><a class="code" href="classCS_1_1Math_1_1Ode45.html" title="Embedded Runge-Kutta 4/5th order ODE solver for non-stiff ODEs.">Ode45</a>
<a name="l00047"></a>00047 {
<a name="l00048"></a>00048 <span class="keyword">public</span>:
<a name="l00049"></a>00049
<a name="l00061"></a>00061 <span class="keyword">template</span><<span class="keyword">typename</span> FuncType, <span class="keyword">typename</span> ArgType>
<a name="l00062"></a><a class="code" href="classCS_1_1Math_1_1Ode45.html#eb07d672fb71fe274e22b6ac768fb30c">00062</a> <span class="keyword">static</span> ArgType <a class="code" href="classCS_1_1Math_1_1Ode45.html#eb07d672fb71fe274e22b6ac768fb30c" title="Step system a single step with step length h.">Step</a> (FuncType& f, ArgType h, ArgType t0, ArgType* y0,
<a name="l00063"></a>00063 ArgType* yout, <span class="keywordtype">size_t</span> size)
<a name="l00064"></a>00064 {
<a name="l00065"></a>00065 <span class="comment">// We need k1-k6</span>
<a name="l00066"></a>00066 <a class="code" href="cssysdef_8h.html#8185c7b21a22716685e7608a1f31a27f" title="Dynamic stack memory allocation.">CS_ALLOC_STACK_ARRAY</a>(ArgType, k1, size);
<a name="l00067"></a>00067 <a class="code" href="cssysdef_8h.html#8185c7b21a22716685e7608a1f31a27f" title="Dynamic stack memory allocation.">CS_ALLOC_STACK_ARRAY</a>(ArgType, k2, size);
<a name="l00068"></a>00068 <a class="code" href="cssysdef_8h.html#8185c7b21a22716685e7608a1f31a27f" title="Dynamic stack memory allocation.">CS_ALLOC_STACK_ARRAY</a>(ArgType, k3, size);
<a name="l00069"></a>00069 <a class="code" href="cssysdef_8h.html#8185c7b21a22716685e7608a1f31a27f" title="Dynamic stack memory allocation.">CS_ALLOC_STACK_ARRAY</a>(ArgType, k4, size);
<a name="l00070"></a>00070 <a class="code" href="cssysdef_8h.html#8185c7b21a22716685e7608a1f31a27f" title="Dynamic stack memory allocation.">CS_ALLOC_STACK_ARRAY</a>(ArgType, k5, size);
<a name="l00071"></a>00071 <a class="code" href="cssysdef_8h.html#8185c7b21a22716685e7608a1f31a27f" title="Dynamic stack memory allocation.">CS_ALLOC_STACK_ARRAY</a>(ArgType, k6, size);
<a name="l00072"></a>00072 <a class="code" href="cssysdef_8h.html#8185c7b21a22716685e7608a1f31a27f" title="Dynamic stack memory allocation.">CS_ALLOC_STACK_ARRAY</a>(ArgType, tmp, size);
<a name="l00073"></a>00073
<a name="l00074"></a>00074 <span class="comment">// k1</span>
<a name="l00075"></a>00075 f (t0, y0, k1, size);
<a name="l00076"></a>00076
<a name="l00077"></a>00077 <span class="comment">// prepare for k2</span>
<a name="l00078"></a>00078 <span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> i = 0; i < size; ++i)
<a name="l00079"></a>00079 {
<a name="l00080"></a>00080 k1[i] *= h;
<a name="l00081"></a>00081 tmp[i] = y0[i] + 0.25*k1[i];
<a name="l00082"></a>00082 }
<a name="l00083"></a>00083
<a name="l00084"></a>00084 <span class="comment">// k2</span>
<a name="l00085"></a>00085 f (t0 + 0.25f*h, tmp, k2, size);
<a name="l00086"></a>00086
<a name="l00087"></a>00087 <span class="comment">// prepare for k3</span>
<a name="l00088"></a>00088 <span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> i = 0; i < size; ++i)
<a name="l00089"></a>00089 {
<a name="l00090"></a>00090 k2[i] *= h;
<a name="l00091"></a>00091 tmp[i] = y0[i] + (3.0/32)*k1[i]
<a name="l00092"></a>00092 + (9.0/32)*k2[i];
<a name="l00093"></a>00093 }
<a name="l00094"></a>00094
<a name="l00095"></a>00095 <span class="comment">// k3</span>
<a name="l00096"></a>00096 f (t0 + (3.0f/8)*h, tmp, k3, size);
<a name="l00097"></a>00097
<a name="l00098"></a>00098 <span class="comment">// prepare for k4</span>
<a name="l00099"></a>00099 <span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> i = 0; i < size; ++i)
<a name="l00100"></a>00100 {
<a name="l00101"></a>00101 k3[i] *= h;
<a name="l00102"></a>00102 tmp[i] = y0[i] + (1932.0/2197)*k1[i]
<a name="l00103"></a>00103 - (7200.0/2197)*k2[i]
<a name="l00104"></a>00104 + (7296.0/2197)*k3[i];
<a name="l00105"></a>00105 }
<a name="l00106"></a>00106
<a name="l00107"></a>00107 <span class="comment">// k4</span>
<a name="l00108"></a>00108 f (t0 + (12.0f/13)*h, tmp, k4, size);
<a name="l00109"></a>00109
<a name="l00110"></a>00110 <span class="comment">// prepare for k5</span>
<a name="l00111"></a>00111 <span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> i = 0; i < size; ++i)
<a name="l00112"></a>00112 {
<a name="l00113"></a>00113 k4[i] *= h;
<a name="l00114"></a>00114 tmp[i] = y0[i] + (439.0/216)*k1[i]
<a name="l00115"></a>00115 - (8.0)*k2[i]
<a name="l00116"></a>00116 + (3680.0/513)*k3[i]
<a name="l00117"></a>00117 - (845.0/4104)*k4[i];
<a name="l00118"></a>00118 }
<a name="l00119"></a>00119
<a name="l00120"></a>00120 <span class="comment">// k5</span>
<a name="l00121"></a>00121 f (t0 + h, tmp, k5, size);
<a name="l00122"></a>00122
<a name="l00123"></a>00123 <span class="comment">// prepare for k6</span>
<a name="l00124"></a>00124 <span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> i = 0; i < size; ++i)
<a name="l00125"></a>00125 {
<a name="l00126"></a>00126 k5[i] *= h;
<a name="l00127"></a>00127 tmp[i] = y0[i] - (8.0/27)*k1[i]
<a name="l00128"></a>00128 + (2.0)*k2[i]
<a name="l00129"></a>00129 - (3544.0/2565)*k3[i]
<a name="l00130"></a>00130 + (1859.0/4104)*k4[i]
<a name="l00131"></a>00131 - (11.0/40)*k5[i];
<a name="l00132"></a>00132 }
<a name="l00133"></a>00133
<a name="l00134"></a>00134 <span class="comment">// k6</span>
<a name="l00135"></a>00135 f (t0 + 0.5f*h, tmp, k6, size);
<a name="l00136"></a>00136
<a name="l00137"></a>00137 ArgType errMag = 0;
<a name="l00138"></a>00138 <span class="comment">// Finally calculate 4th and 5th order result, error term and final result</span>
<a name="l00139"></a>00139 <span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> i = 0; i < size; ++i)
<a name="l00140"></a>00140 {
<a name="l00141"></a>00141
<a name="l00142"></a>00142 k6[i] *= h;
<a name="l00143"></a>00143
<a name="l00144"></a>00144 ArgType y4 = y0[i] + (25.0/216)*k1[i]
<a name="l00145"></a>00145 + (1408.0/2565)*k3[i]
<a name="l00146"></a>00146 + (2197.0/4104)*k4[i]
<a name="l00147"></a>00147 - (1.0/5)*k5[i];
<a name="l00148"></a>00148
<a name="l00149"></a>00149 ArgType y5 = y0[i] + (16.0/315)*k1[i]
<a name="l00150"></a>00150 + (6656.0/12825)*k3[i]
<a name="l00151"></a>00151 + (28561.0/56430)*k4[i]
<a name="l00152"></a>00152 - (9.0f/50)*k5[i]
<a name="l00153"></a>00153 + (2.0/55)*k6[i];
<a name="l00154"></a>00154
<a name="l00155"></a>00155 ArgType yErr = y4 - y5;
<a name="l00156"></a>00156
<a name="l00157"></a>00157 yout[i] = y5 + yErr;
<a name="l00158"></a>00158
<a name="l00159"></a>00159 errMag += yErr*yErr;
<a name="l00160"></a>00160 }
<a name="l00161"></a>00161
<a name="l00162"></a>00162 <span class="keywordflow">return</span> sqrtf (errMag);
<a name="l00163"></a>00163 }
<a name="l00164"></a>00164
<a name="l00165"></a>00165
<a name="l00176"></a>00176 <span class="keyword">template</span><<span class="keyword">typename</span> FuncType, <span class="keyword">typename</span> ArgType>
<a name="l00177"></a><a class="code" href="classCS_1_1Math_1_1Ode45.html#0a9fc38bbee91c52e0edf079dd51acb4">00177</a> <span class="keyword">static</span> <span class="keywordtype">float</span> <a class="code" href="classCS_1_1Math_1_1Ode45.html#eb07d672fb71fe274e22b6ac768fb30c" title="Step system a single step with step length h.">Step</a> (FuncType& f, ArgType h, ArgType t0, <a class="code" href="classcsVector3.html" title="A 3D vector.">csVector3</a> y0,
<a name="l00178"></a>00178 <a class="code" href="classcsVector3.html" title="A 3D vector.">csVector3</a>& yout)
<a name="l00179"></a>00179 {
<a name="l00180"></a>00180 <span class="comment">// We need k1-k6</span>
<a name="l00181"></a>00181 <a class="code" href="classcsVector3.html" title="A 3D vector.">csVector3</a> k1, k2, k3, k4, k5, k6;
<a name="l00182"></a>00182
<a name="l00183"></a>00183 <span class="comment">// k1</span>
<a name="l00184"></a>00184 k1 = h * f (t0, y0);
<a name="l00185"></a>00185
<a name="l00186"></a>00186 <span class="comment">// k2</span>
<a name="l00187"></a>00187 k2 = h * f (t0 + 0.25f*h, y0 + 0.25f*k1);
<a name="l00188"></a>00188
<a name="l00189"></a>00189 <span class="comment">// k3</span>
<a name="l00190"></a>00190 k3 = h * f (t0 + (3.0f/8)*h, y0 + (3.0f/32)*k1
<a name="l00191"></a>00191 + (9.0f/32)*k2);
<a name="l00192"></a>00192
<a name="l00193"></a>00193 <span class="comment">// k4</span>
<a name="l00194"></a>00194 k4 = h * f (t0 + (12.0f/13)*h, y0 + (1932.0f/2197)*k1
<a name="l00195"></a>00195 - (7200.0f/2197)*k2
<a name="l00196"></a>00196 + (7296.0f/2197)*k3);
<a name="l00197"></a>00197
<a name="l00198"></a>00198 <span class="comment">// k5</span>
<a name="l00199"></a>00199 k5 = h * f (t0 + h, y0 + (439.0f/216)*k1
<a name="l00200"></a>00200 - 8.0f*k2
<a name="l00201"></a>00201 + (3680.0f/513)*k3
<a name="l00202"></a>00202 - (845.0f/4104)*k4);
<a name="l00203"></a>00203
<a name="l00204"></a>00204 <span class="comment">// k6</span>
<a name="l00205"></a>00205 k6 = h * f (t0 + 0.5f*h, y0 - (8.0f/27)*k1
<a name="l00206"></a>00206 + (2.0f)*k2
<a name="l00207"></a>00207 - (3544.0f/2565)*k3
<a name="l00208"></a>00208 + (1859.0f/4104)*k4
<a name="l00209"></a>00209 - (11.0f/40)*k5);
<a name="l00210"></a>00210
<a name="l00211"></a>00211 <span class="comment">// Finally calculate 4th and 5th order result, error term and final result</span>
<a name="l00212"></a>00212 <a class="code" href="classcsVector3.html" title="A 3D vector.">csVector3</a> y4 = y0 + (25.0f/216)*k1
<a name="l00213"></a>00213 + (1408.0f/2565)*k3
<a name="l00214"></a>00214 + (2197.0f/4104)*k4
<a name="l00215"></a>00215 - (1.0f/5)*k5;
<a name="l00216"></a>00216
<a name="l00217"></a>00217 <a class="code" href="classcsVector3.html" title="A 3D vector.">csVector3</a> y5 = y0 + (16.0f/315)*k1
<a name="l00218"></a>00218 + (6656.0f/12825)*k3
<a name="l00219"></a>00219 + (28561.0f/56430)*k4
<a name="l00220"></a>00220 - (9.0f/50)*k5
<a name="l00221"></a>00221 + (2.0f/55)*k6;
<a name="l00222"></a>00222
<a name="l00223"></a>00223 <a class="code" href="classcsVector3.html" title="A 3D vector.">csVector3</a> yErr = y4 - y5;
<a name="l00224"></a>00224
<a name="l00225"></a>00225 yout = y5 + yErr;
<a name="l00226"></a>00226
<a name="l00227"></a>00227 <span class="keywordflow">return</span> yErr.<a class="code" href="classcsVector3.html#42b1546a61c4fd76ec17c50592fcad1e" title="Returns the norm of this vector.">Norm</a> ();
<a name="l00228"></a>00228 }
<a name="l00229"></a>00229
<a name="l00240"></a>00240 <span class="keyword">template</span><<span class="keyword">typename</span> FuncType, <span class="keyword">typename</span> ArgType>
<a name="l00241"></a><a class="code" href="classCS_1_1Math_1_1Ode45.html#c8085461e3e8b411a67a99c070cb33b7">00241</a> <span class="keyword">static</span> ArgType <a class="code" href="classCS_1_1Math_1_1Ode45.html#eb07d672fb71fe274e22b6ac768fb30c" title="Step system a single step with step length h.">Step</a> (FuncType& f, ArgType h, ArgType t0, ArgType y0,
<a name="l00242"></a>00242 ArgType& yout)
<a name="l00243"></a>00243 {
<a name="l00244"></a>00244 <span class="comment">// We need k1-k6</span>
<a name="l00245"></a>00245 ArgType k1, k2, k3, k4, k5, k6;
<a name="l00246"></a>00246
<a name="l00247"></a>00247 <span class="comment">// k1</span>
<a name="l00248"></a>00248 k1 = h * f (t0, y0);
<a name="l00249"></a>00249
<a name="l00250"></a>00250 <span class="comment">// k2</span>
<a name="l00251"></a>00251 k2 = h * f (t0 + 0.25f*h, y0 + 0.25f*k1);
<a name="l00252"></a>00252
<a name="l00253"></a>00253 <span class="comment">// k3</span>
<a name="l00254"></a>00254 k3 = h * f (t0 + (3.0f/8)*h, y0 + (3.0f/32)*k1
<a name="l00255"></a>00255 + (9.0f/32)*k2);
<a name="l00256"></a>00256
<a name="l00257"></a>00257 <span class="comment">// k4</span>
<a name="l00258"></a>00258 k4 = h * f (t0 + (12.0f/13)*h, y0 + (1932.0f/2197)*k1
<a name="l00259"></a>00259 - (7200.0f/2197)*k2
<a name="l00260"></a>00260 + (7296.0f/2197)*k3);
<a name="l00261"></a>00261
<a name="l00262"></a>00262 <span class="comment">// k5</span>
<a name="l00263"></a>00263 k5 = h * f (t0 + h, y0 + (439.0f/216)*k1
<a name="l00264"></a>00264 - 8.0f*k2
<a name="l00265"></a>00265 + (3680.0f/513)*k3
<a name="l00266"></a>00266 - (845.0f/4104)*k4);
<a name="l00267"></a>00267
<a name="l00268"></a>00268 <span class="comment">// k6</span>
<a name="l00269"></a>00269 k6 = h * f (t0 + 0.5f*h, y0 - (8.0f/27)*k1
<a name="l00270"></a>00270 + (2.0f)*k2
<a name="l00271"></a>00271 - (3544.0f/2565)*k3
<a name="l00272"></a>00272 + (1859.0f/4104)*k4
<a name="l00273"></a>00273 - (11.0f/40)*k5);
<a name="l00274"></a>00274
<a name="l00275"></a>00275 <span class="comment">// Finally calculate 4th and 5th order result, error term and final result</span>
<a name="l00276"></a>00276 ArgType y4 = y0 + (25.0f/216)*k1
<a name="l00277"></a>00277 + (1408.0f/2565)*k3
<a name="l00278"></a>00278 + (2197.0f/4104)*k4
<a name="l00279"></a>00279 - (1.0f/5)*k5;
<a name="l00280"></a>00280
<a name="l00281"></a>00281 ArgType y5 = y0 + (16.0f/315)*k1
<a name="l00282"></a>00282 + (6656.0f/12825)*k3
<a name="l00283"></a>00283 + (28561.0f/56430)*k4
<a name="l00284"></a>00284 - (9.0f/50)*k5
<a name="l00285"></a>00285 + (2.0f/55)*k6;
<a name="l00286"></a>00286
<a name="l00287"></a>00287 ArgType yErr = y4 - y5;
<a name="l00288"></a>00288
<a name="l00289"></a>00289 yout = y5 + yErr;
<a name="l00290"></a>00290
<a name="l00291"></a>00291 <span class="keywordflow">return</span> yErr;
<a name="l00292"></a>00292 }
<a name="l00293"></a>00293
<a name="l00294"></a>00294 };
<a name="l00295"></a>00295
<a name="l00296"></a>00296
<a name="l00297"></a>00297 }
<a name="l00298"></a>00298 }
<a name="l00299"></a>00299
<a name="l00300"></a>00300 <span class="preprocessor">#endif</span>
</pre></div><hr size="1"><address><small>Generated for Crystal Space 1.2.1 by
<a href="http://www.doxygen.org/index.html">doxygen</a> 1.5.3
</small></address> </div></body> </html>
distrib
>
Mandriva
>
2010.0
>
i586
>
media
>
contrib-release
>
by-pkgid
>
bad97183153701b09df5fae1052b1c30
>
files
>
2095
