Mercurial > hg > de.mpg.mpiwg.itgroup.digilib.plugin

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.html @ 13:cbf34dd4d7e6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
commons-math-2.1 added
author dwinter
date Tue, 04 Jan 2011 10:02:07 +0100
parents
children
comparison
equal deleted inserted replaced
12:970d26a94fb7 13:cbf34dd4d7e6
1 <HTML>
2 <BODY BGCOLOR="white">
3 <PRE>
4 <FONT color="green">001</FONT> /*<a name="line.1"></a>
5 <FONT color="green">002</FONT> * Licensed to the Apache Software Foundation (ASF) under one or more<a name="line.2"></a>
6 <FONT color="green">003</FONT> * contributor license agreements. See the NOTICE file distributed with<a name="line.3"></a>
7 <FONT color="green">004</FONT> * this work for additional information regarding copyright ownership.<a name="line.4"></a>
8 <FONT color="green">005</FONT> * The ASF licenses this file to You under the Apache License, Version 2.0<a name="line.5"></a>
9 <FONT color="green">006</FONT> * (the "License"); you may not use this file except in compliance with<a name="line.6"></a>
10 <FONT color="green">007</FONT> * the License. You may obtain a copy of the License at<a name="line.7"></a>
11 <FONT color="green">008</FONT> *<a name="line.8"></a>
12 <FONT color="green">009</FONT> * http://www.apache.org/licenses/LICENSE-2.0<a name="line.9"></a>
13 <FONT color="green">010</FONT> *<a name="line.10"></a>
14 <FONT color="green">011</FONT> * Unless required by applicable law or agreed to in writing, software<a name="line.11"></a>
15 <FONT color="green">012</FONT> * distributed under the License is distributed on an "AS IS" BASIS,<a name="line.12"></a>
16 <FONT color="green">013</FONT> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.<a name="line.13"></a>
17 <FONT color="green">014</FONT> * See the License for the specific language governing permissions and<a name="line.14"></a>
18 <FONT color="green">015</FONT> * limitations under the License.<a name="line.15"></a>
19 <FONT color="green">016</FONT> */<a name="line.16"></a>
20 <FONT color="green">017</FONT> <a name="line.17"></a>
21 <FONT color="green">018</FONT> package org.apache.commons.math.ode.nonstiff;<a name="line.18"></a>
22 <FONT color="green">019</FONT> <a name="line.19"></a>
23 <FONT color="green">020</FONT> import org.apache.commons.math.ode.AbstractIntegrator;<a name="line.20"></a>
24 <FONT color="green">021</FONT> import org.apache.commons.math.ode.DerivativeException;<a name="line.21"></a>
25 <FONT color="green">022</FONT> import org.apache.commons.math.ode.FirstOrderDifferentialEquations;<a name="line.22"></a>
26 <FONT color="green">023</FONT> import org.apache.commons.math.ode.IntegratorException;<a name="line.23"></a>
27 <FONT color="green">024</FONT> <a name="line.24"></a>
28 <FONT color="green">025</FONT> /**<a name="line.25"></a>
29 <FONT color="green">026</FONT> * This abstract class holds the common part of all adaptive<a name="line.26"></a>
30 <FONT color="green">027</FONT> * stepsize integrators for Ordinary Differential Equations.<a name="line.27"></a>
31 <FONT color="green">028</FONT> *<a name="line.28"></a>
32 <FONT color="green">029</FONT> * &lt;p&gt;These algorithms perform integration with stepsize control, which<a name="line.29"></a>
33 <FONT color="green">030</FONT> * means the user does not specify the integration step but rather a<a name="line.30"></a>
34 <FONT color="green">031</FONT> * tolerance on error. The error threshold is computed as<a name="line.31"></a>
35 <FONT color="green">032</FONT> * &lt;pre&gt;<a name="line.32"></a>
36 <FONT color="green">033</FONT> * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))<a name="line.33"></a>
37 <FONT color="green">034</FONT> * &lt;/pre&gt;<a name="line.34"></a>
38 <FONT color="green">035</FONT> * where absTol_i is the absolute tolerance for component i of the<a name="line.35"></a>
39 <FONT color="green">036</FONT> * state vector and relTol_i is the relative tolerance for the same<a name="line.36"></a>
40 <FONT color="green">037</FONT> * component. The user can also use only two scalar values absTol and<a name="line.37"></a>
41 <FONT color="green">038</FONT> * relTol which will be used for all components.&lt;/p&gt;<a name="line.38"></a>
42 <FONT color="green">039</FONT> *<a name="line.39"></a>
43 <FONT color="green">040</FONT> * &lt;p&gt;If the estimated error for ym+1 is such that<a name="line.40"></a>
44 <FONT color="green">041</FONT> * &lt;pre&gt;<a name="line.41"></a>
45 <FONT color="green">042</FONT> * sqrt((sum (errEst_i / threshold_i)^2 ) / n) &lt; 1<a name="line.42"></a>
46 <FONT color="green">043</FONT> * &lt;/pre&gt;<a name="line.43"></a>
47 <FONT color="green">044</FONT> *<a name="line.44"></a>
48 <FONT color="green">045</FONT> * (where n is the state vector dimension) then the step is accepted,<a name="line.45"></a>
49 <FONT color="green">046</FONT> * otherwise the step is rejected and a new attempt is made with a new<a name="line.46"></a>
50 <FONT color="green">047</FONT> * stepsize.&lt;/p&gt;<a name="line.47"></a>
51 <FONT color="green">048</FONT> *<a name="line.48"></a>
52 <FONT color="green">049</FONT> * @version $Revision: 811827 $ $Date: 2009-09-06 11:32:50 -0400 (Sun, 06 Sep 2009) $<a name="line.49"></a>
53 <FONT color="green">050</FONT> * @since 1.2<a name="line.50"></a>
54 <FONT color="green">051</FONT> *<a name="line.51"></a>
55 <FONT color="green">052</FONT> */<a name="line.52"></a>
56 <FONT color="green">053</FONT> <a name="line.53"></a>
57 <FONT color="green">054</FONT> public abstract class AdaptiveStepsizeIntegrator<a name="line.54"></a>
58 <FONT color="green">055</FONT> extends AbstractIntegrator {<a name="line.55"></a>
59 <FONT color="green">056</FONT> <a name="line.56"></a>
60 <FONT color="green">057</FONT> /** Allowed absolute scalar error. */<a name="line.57"></a>
61 <FONT color="green">058</FONT> protected final double scalAbsoluteTolerance;<a name="line.58"></a>
62 <FONT color="green">059</FONT> <a name="line.59"></a>
63 <FONT color="green">060</FONT> /** Allowed relative scalar error. */<a name="line.60"></a>
64 <FONT color="green">061</FONT> protected final double scalRelativeTolerance;<a name="line.61"></a>
65 <FONT color="green">062</FONT> <a name="line.62"></a>
66 <FONT color="green">063</FONT> /** Allowed absolute vectorial error. */<a name="line.63"></a>
67 <FONT color="green">064</FONT> protected final double[] vecAbsoluteTolerance;<a name="line.64"></a>
68 <FONT color="green">065</FONT> <a name="line.65"></a>
69 <FONT color="green">066</FONT> /** Allowed relative vectorial error. */<a name="line.66"></a>
70 <FONT color="green">067</FONT> protected final double[] vecRelativeTolerance;<a name="line.67"></a>
71 <FONT color="green">068</FONT> <a name="line.68"></a>
72 <FONT color="green">069</FONT> /** User supplied initial step. */<a name="line.69"></a>
73 <FONT color="green">070</FONT> private double initialStep;<a name="line.70"></a>
74 <FONT color="green">071</FONT> <a name="line.71"></a>
75 <FONT color="green">072</FONT> /** Minimal step. */<a name="line.72"></a>
76 <FONT color="green">073</FONT> private final double minStep;<a name="line.73"></a>
77 <FONT color="green">074</FONT> <a name="line.74"></a>
78 <FONT color="green">075</FONT> /** Maximal step. */<a name="line.75"></a>
79 <FONT color="green">076</FONT> private final double maxStep;<a name="line.76"></a>
80 <FONT color="green">077</FONT> <a name="line.77"></a>
81 <FONT color="green">078</FONT> /** Build an integrator with the given stepsize bounds.<a name="line.78"></a>
82 <FONT color="green">079</FONT> * The default step handler does nothing.<a name="line.79"></a>
83 <FONT color="green">080</FONT> * @param name name of the method<a name="line.80"></a>
84 <FONT color="green">081</FONT> * @param minStep minimal step (must be positive even for backward<a name="line.81"></a>
85 <FONT color="green">082</FONT> * integration), the last step can be smaller than this<a name="line.82"></a>
86 <FONT color="green">083</FONT> * @param maxStep maximal step (must be positive even for backward<a name="line.83"></a>
87 <FONT color="green">084</FONT> * integration)<a name="line.84"></a>
88 <FONT color="green">085</FONT> * @param scalAbsoluteTolerance allowed absolute error<a name="line.85"></a>
89 <FONT color="green">086</FONT> * @param scalRelativeTolerance allowed relative error<a name="line.86"></a>
90 <FONT color="green">087</FONT> */<a name="line.87"></a>
91 <FONT color="green">088</FONT> public AdaptiveStepsizeIntegrator(final String name,<a name="line.88"></a>
92 <FONT color="green">089</FONT> final double minStep, final double maxStep,<a name="line.89"></a>
93 <FONT color="green">090</FONT> final double scalAbsoluteTolerance,<a name="line.90"></a>
94 <FONT color="green">091</FONT> final double scalRelativeTolerance) {<a name="line.91"></a>
95 <FONT color="green">092</FONT> <a name="line.92"></a>
96 <FONT color="green">093</FONT> super(name);<a name="line.93"></a>
97 <FONT color="green">094</FONT> <a name="line.94"></a>
98 <FONT color="green">095</FONT> this.minStep = Math.abs(minStep);<a name="line.95"></a>
99 <FONT color="green">096</FONT> this.maxStep = Math.abs(maxStep);<a name="line.96"></a>
100 <FONT color="green">097</FONT> this.initialStep = -1.0;<a name="line.97"></a>
101 <FONT color="green">098</FONT> <a name="line.98"></a>
102 <FONT color="green">099</FONT> this.scalAbsoluteTolerance = scalAbsoluteTolerance;<a name="line.99"></a>
103 <FONT color="green">100</FONT> this.scalRelativeTolerance = scalRelativeTolerance;<a name="line.100"></a>
104 <FONT color="green">101</FONT> this.vecAbsoluteTolerance = null;<a name="line.101"></a>
105 <FONT color="green">102</FONT> this.vecRelativeTolerance = null;<a name="line.102"></a>
106 <FONT color="green">103</FONT> <a name="line.103"></a>
107 <FONT color="green">104</FONT> resetInternalState();<a name="line.104"></a>
108 <FONT color="green">105</FONT> <a name="line.105"></a>
109 <FONT color="green">106</FONT> }<a name="line.106"></a>
110 <FONT color="green">107</FONT> <a name="line.107"></a>
111 <FONT color="green">108</FONT> /** Build an integrator with the given stepsize bounds.<a name="line.108"></a>
112 <FONT color="green">109</FONT> * The default step handler does nothing.<a name="line.109"></a>
113 <FONT color="green">110</FONT> * @param name name of the method<a name="line.110"></a>
114 <FONT color="green">111</FONT> * @param minStep minimal step (must be positive even for backward<a name="line.111"></a>
115 <FONT color="green">112</FONT> * integration), the last step can be smaller than this<a name="line.112"></a>
116 <FONT color="green">113</FONT> * @param maxStep maximal step (must be positive even for backward<a name="line.113"></a>
117 <FONT color="green">114</FONT> * integration)<a name="line.114"></a>
118 <FONT color="green">115</FONT> * @param vecAbsoluteTolerance allowed absolute error<a name="line.115"></a>
119 <FONT color="green">116</FONT> * @param vecRelativeTolerance allowed relative error<a name="line.116"></a>
120 <FONT color="green">117</FONT> */<a name="line.117"></a>
121 <FONT color="green">118</FONT> public AdaptiveStepsizeIntegrator(final String name,<a name="line.118"></a>
122 <FONT color="green">119</FONT> final double minStep, final double maxStep,<a name="line.119"></a>
123 <FONT color="green">120</FONT> final double[] vecAbsoluteTolerance,<a name="line.120"></a>
124 <FONT color="green">121</FONT> final double[] vecRelativeTolerance) {<a name="line.121"></a>
125 <FONT color="green">122</FONT> <a name="line.122"></a>
126 <FONT color="green">123</FONT> super(name);<a name="line.123"></a>
127 <FONT color="green">124</FONT> <a name="line.124"></a>
128 <FONT color="green">125</FONT> this.minStep = minStep;<a name="line.125"></a>
129 <FONT color="green">126</FONT> this.maxStep = maxStep;<a name="line.126"></a>
130 <FONT color="green">127</FONT> this.initialStep = -1.0;<a name="line.127"></a>
131 <FONT color="green">128</FONT> <a name="line.128"></a>
132 <FONT color="green">129</FONT> this.scalAbsoluteTolerance = 0;<a name="line.129"></a>
133 <FONT color="green">130</FONT> this.scalRelativeTolerance = 0;<a name="line.130"></a>
134 <FONT color="green">131</FONT> this.vecAbsoluteTolerance = vecAbsoluteTolerance.clone();<a name="line.131"></a>
135 <FONT color="green">132</FONT> this.vecRelativeTolerance = vecRelativeTolerance.clone();<a name="line.132"></a>
136 <FONT color="green">133</FONT> <a name="line.133"></a>
137 <FONT color="green">134</FONT> resetInternalState();<a name="line.134"></a>
138 <FONT color="green">135</FONT> <a name="line.135"></a>
139 <FONT color="green">136</FONT> }<a name="line.136"></a>
140 <FONT color="green">137</FONT> <a name="line.137"></a>
141 <FONT color="green">138</FONT> /** Set the initial step size.<a name="line.138"></a>
142 <FONT color="green">139</FONT> * &lt;p&gt;This method allows the user to specify an initial positive<a name="line.139"></a>
143 <FONT color="green">140</FONT> * step size instead of letting the integrator guess it by<a name="line.140"></a>
144 <FONT color="green">141</FONT> * itself. If this method is not called before integration is<a name="line.141"></a>
145 <FONT color="green">142</FONT> * started, the initial step size will be estimated by the<a name="line.142"></a>
146 <FONT color="green">143</FONT> * integrator.&lt;/p&gt;<a name="line.143"></a>
147 <FONT color="green">144</FONT> * @param initialStepSize initial step size to use (must be positive even<a name="line.144"></a>
148 <FONT color="green">145</FONT> * for backward integration ; providing a negative value or a value<a name="line.145"></a>
149 <FONT color="green">146</FONT> * outside of the min/max step interval will lead the integrator to<a name="line.146"></a>
150 <FONT color="green">147</FONT> * ignore the value and compute the initial step size by itself)<a name="line.147"></a>
151 <FONT color="green">148</FONT> */<a name="line.148"></a>
152 <FONT color="green">149</FONT> public void setInitialStepSize(final double initialStepSize) {<a name="line.149"></a>
153 <FONT color="green">150</FONT> if ((initialStepSize &lt; minStep) || (initialStepSize &gt; maxStep)) {<a name="line.150"></a>
154 <FONT color="green">151</FONT> initialStep = -1.0;<a name="line.151"></a>
155 <FONT color="green">152</FONT> } else {<a name="line.152"></a>
156 <FONT color="green">153</FONT> initialStep = initialStepSize;<a name="line.153"></a>
157 <FONT color="green">154</FONT> }<a name="line.154"></a>
158 <FONT color="green">155</FONT> }<a name="line.155"></a>
159 <FONT color="green">156</FONT> <a name="line.156"></a>
160 <FONT color="green">157</FONT> /** Perform some sanity checks on the integration parameters.<a name="line.157"></a>
161 <FONT color="green">158</FONT> * @param equations differential equations set<a name="line.158"></a>
162 <FONT color="green">159</FONT> * @param t0 start time<a name="line.159"></a>
163 <FONT color="green">160</FONT> * @param y0 state vector at t0<a name="line.160"></a>
164 <FONT color="green">161</FONT> * @param t target time for the integration<a name="line.161"></a>
165 <FONT color="green">162</FONT> * @param y placeholder where to put the state vector<a name="line.162"></a>
166 <FONT color="green">163</FONT> * @exception IntegratorException if some inconsistency is detected<a name="line.163"></a>
167 <FONT color="green">164</FONT> */<a name="line.164"></a>
168 <FONT color="green">165</FONT> @Override<a name="line.165"></a>
169 <FONT color="green">166</FONT> protected void sanityChecks(final FirstOrderDifferentialEquations equations,<a name="line.166"></a>
170 <FONT color="green">167</FONT> final double t0, final double[] y0,<a name="line.167"></a>
171 <FONT color="green">168</FONT> final double t, final double[] y)<a name="line.168"></a>
172 <FONT color="green">169</FONT> throws IntegratorException {<a name="line.169"></a>
173 <FONT color="green">170</FONT> <a name="line.170"></a>
174 <FONT color="green">171</FONT> super.sanityChecks(equations, t0, y0, t, y);<a name="line.171"></a>
175 <FONT color="green">172</FONT> <a name="line.172"></a>
176 <FONT color="green">173</FONT> if ((vecAbsoluteTolerance != null) &amp;&amp; (vecAbsoluteTolerance.length != y0.length)) {<a name="line.173"></a>
177 <FONT color="green">174</FONT> throw new IntegratorException(<a name="line.174"></a>
178 <FONT color="green">175</FONT> "dimensions mismatch: state vector has dimension {0}," +<a name="line.175"></a>
179 <FONT color="green">176</FONT> " absolute tolerance vector has dimension {1}",<a name="line.176"></a>
180 <FONT color="green">177</FONT> y0.length, vecAbsoluteTolerance.length);<a name="line.177"></a>
181 <FONT color="green">178</FONT> }<a name="line.178"></a>
182 <FONT color="green">179</FONT> <a name="line.179"></a>
183 <FONT color="green">180</FONT> if ((vecRelativeTolerance != null) &amp;&amp; (vecRelativeTolerance.length != y0.length)) {<a name="line.180"></a>
184 <FONT color="green">181</FONT> throw new IntegratorException(<a name="line.181"></a>
185 <FONT color="green">182</FONT> "dimensions mismatch: state vector has dimension {0}," +<a name="line.182"></a>
186 <FONT color="green">183</FONT> " relative tolerance vector has dimension {1}",<a name="line.183"></a>
187 <FONT color="green">184</FONT> y0.length, vecRelativeTolerance.length);<a name="line.184"></a>
188 <FONT color="green">185</FONT> }<a name="line.185"></a>
189 <FONT color="green">186</FONT> <a name="line.186"></a>
190 <FONT color="green">187</FONT> }<a name="line.187"></a>
191 <FONT color="green">188</FONT> <a name="line.188"></a>
192 <FONT color="green">189</FONT> /** Initialize the integration step.<a name="line.189"></a>
193 <FONT color="green">190</FONT> * @param equations differential equations set<a name="line.190"></a>
194 <FONT color="green">191</FONT> * @param forward forward integration indicator<a name="line.191"></a>
195 <FONT color="green">192</FONT> * @param order order of the method<a name="line.192"></a>
196 <FONT color="green">193</FONT> * @param scale scaling vector for the state vector<a name="line.193"></a>
197 <FONT color="green">194</FONT> * @param t0 start time<a name="line.194"></a>
198 <FONT color="green">195</FONT> * @param y0 state vector at t0<a name="line.195"></a>
199 <FONT color="green">196</FONT> * @param yDot0 first time derivative of y0<a name="line.196"></a>
200 <FONT color="green">197</FONT> * @param y1 work array for a state vector<a name="line.197"></a>
201 <FONT color="green">198</FONT> * @param yDot1 work array for the first time derivative of y1<a name="line.198"></a>
202 <FONT color="green">199</FONT> * @return first integration step<a name="line.199"></a>
203 <FONT color="green">200</FONT> * @exception DerivativeException this exception is propagated to<a name="line.200"></a>
204 <FONT color="green">201</FONT> * the caller if the underlying user function triggers one<a name="line.201"></a>
205 <FONT color="green">202</FONT> */<a name="line.202"></a>
206 <FONT color="green">203</FONT> public double initializeStep(final FirstOrderDifferentialEquations equations,<a name="line.203"></a>
207 <FONT color="green">204</FONT> final boolean forward, final int order, final double[] scale,<a name="line.204"></a>
208 <FONT color="green">205</FONT> final double t0, final double[] y0, final double[] yDot0,<a name="line.205"></a>
209 <FONT color="green">206</FONT> final double[] y1, final double[] yDot1)<a name="line.206"></a>
210 <FONT color="green">207</FONT> throws DerivativeException {<a name="line.207"></a>
211 <FONT color="green">208</FONT> <a name="line.208"></a>
212 <FONT color="green">209</FONT> if (initialStep &gt; 0) {<a name="line.209"></a>
213 <FONT color="green">210</FONT> // use the user provided value<a name="line.210"></a>
214 <FONT color="green">211</FONT> return forward ? initialStep : -initialStep;<a name="line.211"></a>
215 <FONT color="green">212</FONT> }<a name="line.212"></a>
216 <FONT color="green">213</FONT> <a name="line.213"></a>
217 <FONT color="green">214</FONT> // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||<a name="line.214"></a>
218 <FONT color="green">215</FONT> // this guess will be used to perform an Euler step<a name="line.215"></a>
219 <FONT color="green">216</FONT> double ratio;<a name="line.216"></a>
220 <FONT color="green">217</FONT> double yOnScale2 = 0;<a name="line.217"></a>
221 <FONT color="green">218</FONT> double yDotOnScale2 = 0;<a name="line.218"></a>
222 <FONT color="green">219</FONT> for (int j = 0; j &lt; y0.length; ++j) {<a name="line.219"></a>
223 <FONT color="green">220</FONT> ratio = y0[j] / scale[j];<a name="line.220"></a>
224 <FONT color="green">221</FONT> yOnScale2 += ratio * ratio;<a name="line.221"></a>
225 <FONT color="green">222</FONT> ratio = yDot0[j] / scale[j];<a name="line.222"></a>
226 <FONT color="green">223</FONT> yDotOnScale2 += ratio * ratio;<a name="line.223"></a>
227 <FONT color="green">224</FONT> }<a name="line.224"></a>
228 <FONT color="green">225</FONT> <a name="line.225"></a>
229 <FONT color="green">226</FONT> double h = ((yOnScale2 &lt; 1.0e-10) || (yDotOnScale2 &lt; 1.0e-10)) ?<a name="line.226"></a>
230 <FONT color="green">227</FONT> 1.0e-6 : (0.01 * Math.sqrt(yOnScale2 / yDotOnScale2));<a name="line.227"></a>
231 <FONT color="green">228</FONT> if (! forward) {<a name="line.228"></a>
232 <FONT color="green">229</FONT> h = -h;<a name="line.229"></a>
233 <FONT color="green">230</FONT> }<a name="line.230"></a>
234 <FONT color="green">231</FONT> <a name="line.231"></a>
235 <FONT color="green">232</FONT> // perform an Euler step using the preceding rough guess<a name="line.232"></a>
236 <FONT color="green">233</FONT> for (int j = 0; j &lt; y0.length; ++j) {<a name="line.233"></a>
237 <FONT color="green">234</FONT> y1[j] = y0[j] + h * yDot0[j];<a name="line.234"></a>
238 <FONT color="green">235</FONT> }<a name="line.235"></a>
239 <FONT color="green">236</FONT> computeDerivatives(t0 + h, y1, yDot1);<a name="line.236"></a>
240 <FONT color="green">237</FONT> <a name="line.237"></a>
241 <FONT color="green">238</FONT> // estimate the second derivative of the solution<a name="line.238"></a>
242 <FONT color="green">239</FONT> double yDDotOnScale = 0;<a name="line.239"></a>
243 <FONT color="green">240</FONT> for (int j = 0; j &lt; y0.length; ++j) {<a name="line.240"></a>
244 <FONT color="green">241</FONT> ratio = (yDot1[j] - yDot0[j]) / scale[j];<a name="line.241"></a>
245 <FONT color="green">242</FONT> yDDotOnScale += ratio * ratio;<a name="line.242"></a>
246 <FONT color="green">243</FONT> }<a name="line.243"></a>
247 <FONT color="green">244</FONT> yDDotOnScale = Math.sqrt(yDDotOnScale) / h;<a name="line.244"></a>
248 <FONT color="green">245</FONT> <a name="line.245"></a>
249 <FONT color="green">246</FONT> // step size is computed such that<a name="line.246"></a>
250 <FONT color="green">247</FONT> // h^order * max (||y'/tol||, ||y''/tol||) = 0.01<a name="line.247"></a>
251 <FONT color="green">248</FONT> final double maxInv2 = Math.max(Math.sqrt(yDotOnScale2), yDDotOnScale);<a name="line.248"></a>
252 <FONT color="green">249</FONT> final double h1 = (maxInv2 &lt; 1.0e-15) ?<a name="line.249"></a>
253 <FONT color="green">250</FONT> Math.max(1.0e-6, 0.001 * Math.abs(h)) :<a name="line.250"></a>
254 <FONT color="green">251</FONT> Math.pow(0.01 / maxInv2, 1.0 / order);<a name="line.251"></a>
255 <FONT color="green">252</FONT> h = Math.min(100.0 * Math.abs(h), h1);<a name="line.252"></a>
256 <FONT color="green">253</FONT> h = Math.max(h, 1.0e-12 * Math.abs(t0)); // avoids cancellation when computing t1 - t0<a name="line.253"></a>
257 <FONT color="green">254</FONT> if (h &lt; getMinStep()) {<a name="line.254"></a>
258 <FONT color="green">255</FONT> h = getMinStep();<a name="line.255"></a>
259 <FONT color="green">256</FONT> }<a name="line.256"></a>
260 <FONT color="green">257</FONT> if (h &gt; getMaxStep()) {<a name="line.257"></a>
261 <FONT color="green">258</FONT> h = getMaxStep();<a name="line.258"></a>
262 <FONT color="green">259</FONT> }<a name="line.259"></a>
263 <FONT color="green">260</FONT> if (! forward) {<a name="line.260"></a>
264 <FONT color="green">261</FONT> h = -h;<a name="line.261"></a>
265 <FONT color="green">262</FONT> }<a name="line.262"></a>
266 <FONT color="green">263</FONT> <a name="line.263"></a>
267 <FONT color="green">264</FONT> return h;<a name="line.264"></a>
268 <FONT color="green">265</FONT> <a name="line.265"></a>
269 <FONT color="green">266</FONT> }<a name="line.266"></a>
270 <FONT color="green">267</FONT> <a name="line.267"></a>
271 <FONT color="green">268</FONT> /** Filter the integration step.<a name="line.268"></a>
272 <FONT color="green">269</FONT> * @param h signed step<a name="line.269"></a>
273 <FONT color="green">270</FONT> * @param forward forward integration indicator<a name="line.270"></a>
274 <FONT color="green">271</FONT> * @param acceptSmall if true, steps smaller than the minimal value<a name="line.271"></a>
275 <FONT color="green">272</FONT> * are silently increased up to this value, if false such small<a name="line.272"></a>
276 <FONT color="green">273</FONT> * steps generate an exception<a name="line.273"></a>
277 <FONT color="green">274</FONT> * @return a bounded integration step (h if no bound is reach, or a bounded value)<a name="line.274"></a>
278 <FONT color="green">275</FONT> * @exception IntegratorException if the step is too small and acceptSmall is false<a name="line.275"></a>
279 <FONT color="green">276</FONT> */<a name="line.276"></a>
280 <FONT color="green">277</FONT> protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)<a name="line.277"></a>
281 <FONT color="green">278</FONT> throws IntegratorException {<a name="line.278"></a>
282 <FONT color="green">279</FONT> <a name="line.279"></a>
283 <FONT color="green">280</FONT> double filteredH = h;<a name="line.280"></a>
284 <FONT color="green">281</FONT> if (Math.abs(h) &lt; minStep) {<a name="line.281"></a>
285 <FONT color="green">282</FONT> if (acceptSmall) {<a name="line.282"></a>
286 <FONT color="green">283</FONT> filteredH = forward ? minStep : -minStep;<a name="line.283"></a>
287 <FONT color="green">284</FONT> } else {<a name="line.284"></a>
288 <FONT color="green">285</FONT> throw new IntegratorException(<a name="line.285"></a>
289 <FONT color="green">286</FONT> "minimal step size ({0,number,0.00E00}) reached, integration needs {1,number,0.00E00}",<a name="line.286"></a>
290 <FONT color="green">287</FONT> minStep, Math.abs(h));<a name="line.287"></a>
291 <FONT color="green">288</FONT> }<a name="line.288"></a>
292 <FONT color="green">289</FONT> }<a name="line.289"></a>
293 <FONT color="green">290</FONT> <a name="line.290"></a>
294 <FONT color="green">291</FONT> if (filteredH &gt; maxStep) {<a name="line.291"></a>
295 <FONT color="green">292</FONT> filteredH = maxStep;<a name="line.292"></a>
296 <FONT color="green">293</FONT> } else if (filteredH &lt; -maxStep) {<a name="line.293"></a>
297 <FONT color="green">294</FONT> filteredH = -maxStep;<a name="line.294"></a>
298 <FONT color="green">295</FONT> }<a name="line.295"></a>
299 <FONT color="green">296</FONT> <a name="line.296"></a>
300 <FONT color="green">297</FONT> return filteredH;<a name="line.297"></a>
301 <FONT color="green">298</FONT> <a name="line.298"></a>
302 <FONT color="green">299</FONT> }<a name="line.299"></a>
303 <FONT color="green">300</FONT> <a name="line.300"></a>
304 <FONT color="green">301</FONT> /** {@inheritDoc} */<a name="line.301"></a>
305 <FONT color="green">302</FONT> public abstract double integrate (FirstOrderDifferentialEquations equations,<a name="line.302"></a>
306 <FONT color="green">303</FONT> double t0, double[] y0,<a name="line.303"></a>
307 <FONT color="green">304</FONT> double t, double[] y)<a name="line.304"></a>
308 <FONT color="green">305</FONT> throws DerivativeException, IntegratorException;<a name="line.305"></a>
309 <FONT color="green">306</FONT> <a name="line.306"></a>
310 <FONT color="green">307</FONT> /** {@inheritDoc} */<a name="line.307"></a>
311 <FONT color="green">308</FONT> @Override<a name="line.308"></a>
312 <FONT color="green">309</FONT> public double getCurrentStepStart() {<a name="line.309"></a>
313 <FONT color="green">310</FONT> return stepStart;<a name="line.310"></a>
314 <FONT color="green">311</FONT> }<a name="line.311"></a>
315 <FONT color="green">312</FONT> <a name="line.312"></a>
316 <FONT color="green">313</FONT> /** Reset internal state to dummy values. */<a name="line.313"></a>
317 <FONT color="green">314</FONT> protected void resetInternalState() {<a name="line.314"></a>
318 <FONT color="green">315</FONT> stepStart = Double.NaN;<a name="line.315"></a>
319 <FONT color="green">316</FONT> stepSize = Math.sqrt(minStep * maxStep);<a name="line.316"></a>
320 <FONT color="green">317</FONT> }<a name="line.317"></a>
321 <FONT color="green">318</FONT> <a name="line.318"></a>
322 <FONT color="green">319</FONT> /** Get the minimal step.<a name="line.319"></a>
323 <FONT color="green">320</FONT> * @return minimal step<a name="line.320"></a>
324 <FONT color="green">321</FONT> */<a name="line.321"></a>
325 <FONT color="green">322</FONT> public double getMinStep() {<a name="line.322"></a>
326 <FONT color="green">323</FONT> return minStep;<a name="line.323"></a>
327 <FONT color="green">324</FONT> }<a name="line.324"></a>
328 <FONT color="green">325</FONT> <a name="line.325"></a>
329 <FONT color="green">326</FONT> /** Get the maximal step.<a name="line.326"></a>
330 <FONT color="green">327</FONT> * @return maximal step<a name="line.327"></a>
331 <FONT color="green">328</FONT> */<a name="line.328"></a>
332 <FONT color="green">329</FONT> public double getMaxStep() {<a name="line.329"></a>
333 <FONT color="green">330</FONT> return maxStep;<a name="line.330"></a>
334 <FONT color="green">331</FONT> }<a name="line.331"></a>
335 <FONT color="green">332</FONT> <a name="line.332"></a>
336 <FONT color="green">333</FONT> }<a name="line.333"></a>
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397 </PRE>
398 </BODY>
399 </HTML>