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

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/ode/MultistepIntegrator.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;<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.MathRuntimeException;<a name="line.20"></a>
24 <FONT color="green">021</FONT> import org.apache.commons.math.linear.Array2DRowRealMatrix;<a name="line.21"></a>
25 <FONT color="green">022</FONT> import org.apache.commons.math.linear.RealMatrix;<a name="line.22"></a>
26 <FONT color="green">023</FONT> import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;<a name="line.23"></a>
27 <FONT color="green">024</FONT> import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;<a name="line.24"></a>
28 <FONT color="green">025</FONT> import org.apache.commons.math.ode.sampling.StepHandler;<a name="line.25"></a>
29 <FONT color="green">026</FONT> import org.apache.commons.math.ode.sampling.StepInterpolator;<a name="line.26"></a>
30 <FONT color="green">027</FONT> <a name="line.27"></a>
31 <FONT color="green">028</FONT> /**<a name="line.28"></a>
32 <FONT color="green">029</FONT> * This class is the base class for multistep integrators for Ordinary<a name="line.29"></a>
33 <FONT color="green">030</FONT> * Differential Equations.<a name="line.30"></a>
34 <FONT color="green">031</FONT> * &lt;p&gt;We define scaled derivatives s&lt;sub&gt;i&lt;/sub&gt;(n) at step n 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> * s&lt;sub&gt;1&lt;/sub&gt;(n) = h y'&lt;sub&gt;n&lt;/sub&gt; for first derivative<a name="line.33"></a>
37 <FONT color="green">034</FONT> * s&lt;sub&gt;2&lt;/sub&gt;(n) = h&lt;sup&gt;2&lt;/sup&gt;/2 y''&lt;sub&gt;n&lt;/sub&gt; for second derivative<a name="line.34"></a>
38 <FONT color="green">035</FONT> * s&lt;sub&gt;3&lt;/sub&gt;(n) = h&lt;sup&gt;3&lt;/sup&gt;/6 y'''&lt;sub&gt;n&lt;/sub&gt; for third derivative<a name="line.35"></a>
39 <FONT color="green">036</FONT> * ...<a name="line.36"></a>
40 <FONT color="green">037</FONT> * s&lt;sub&gt;k&lt;/sub&gt;(n) = h&lt;sup&gt;k&lt;/sup&gt;/k! y(k)&lt;sub&gt;n&lt;/sub&gt; for k&lt;sup&gt;th&lt;/sup&gt; derivative<a name="line.37"></a>
41 <FONT color="green">038</FONT> * &lt;/pre&gt;&lt;/p&gt;<a name="line.38"></a>
42 <FONT color="green">039</FONT> * &lt;p&gt;Rather than storing several previous steps separately, this implementation uses<a name="line.39"></a>
43 <FONT color="green">040</FONT> * the Nordsieck vector with higher degrees scaled derivatives all taken at the same<a name="line.40"></a>
44 <FONT color="green">041</FONT> * step (y&lt;sub&gt;n&lt;/sub&gt;, s&lt;sub&gt;1&lt;/sub&gt;(n) and r&lt;sub&gt;n&lt;/sub&gt;) where r&lt;sub&gt;n&lt;/sub&gt; is defined as:<a name="line.41"></a>
45 <FONT color="green">042</FONT> * &lt;pre&gt;<a name="line.42"></a>
46 <FONT color="green">043</FONT> * r&lt;sub&gt;n&lt;/sub&gt; = [ s&lt;sub&gt;2&lt;/sub&gt;(n), s&lt;sub&gt;3&lt;/sub&gt;(n) ... s&lt;sub&gt;k&lt;/sub&gt;(n) ]&lt;sup&gt;T&lt;/sup&gt;<a name="line.43"></a>
47 <FONT color="green">044</FONT> * &lt;/pre&gt;<a name="line.44"></a>
48 <FONT color="green">045</FONT> * (we omit the k index in the notation for clarity)&lt;/p&gt;<a name="line.45"></a>
49 <FONT color="green">046</FONT> * &lt;p&gt;<a name="line.46"></a>
50 <FONT color="green">047</FONT> * Multistep integrators with Nordsieck representation are highly sensitive to<a name="line.47"></a>
51 <FONT color="green">048</FONT> * large step changes because when the step is multiplied by a factor a, the<a name="line.48"></a>
52 <FONT color="green">049</FONT> * k&lt;sup&gt;th&lt;/sup&gt; component of the Nordsieck vector is multiplied by a&lt;sup&gt;k&lt;/sup&gt;<a name="line.49"></a>
53 <FONT color="green">050</FONT> * and the last components are the least accurate ones. The default max growth<a name="line.50"></a>
54 <FONT color="green">051</FONT> * factor is therefore set to a quite low value: 2&lt;sup&gt;1/order&lt;/sup&gt;.<a name="line.51"></a>
55 <FONT color="green">052</FONT> * &lt;/p&gt;<a name="line.52"></a>
56 <FONT color="green">053</FONT> *<a name="line.53"></a>
57 <FONT color="green">054</FONT> * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator<a name="line.54"></a>
58 <FONT color="green">055</FONT> * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator<a name="line.55"></a>
59 <FONT color="green">056</FONT> * @version $Revision: 811827 $ $Date: 2009-09-06 11:32:50 -0400 (Sun, 06 Sep 2009) $<a name="line.56"></a>
60 <FONT color="green">057</FONT> * @since 2.0<a name="line.57"></a>
61 <FONT color="green">058</FONT> */<a name="line.58"></a>
62 <FONT color="green">059</FONT> public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {<a name="line.59"></a>
63 <FONT color="green">060</FONT> <a name="line.60"></a>
64 <FONT color="green">061</FONT> /** First scaled derivative (h y'). */<a name="line.61"></a>
65 <FONT color="green">062</FONT> protected double[] scaled;<a name="line.62"></a>
66 <FONT color="green">063</FONT> <a name="line.63"></a>
67 <FONT color="green">064</FONT> /** Nordsieck matrix of the higher scaled derivatives.<a name="line.64"></a>
68 <FONT color="green">065</FONT> * &lt;p&gt;(h&lt;sup&gt;2&lt;/sup&gt;/2 y'', h&lt;sup&gt;3&lt;/sup&gt;/6 y''' ..., h&lt;sup&gt;k&lt;/sup&gt;/k! y(k))&lt;/p&gt;<a name="line.65"></a>
69 <FONT color="green">066</FONT> */<a name="line.66"></a>
70 <FONT color="green">067</FONT> protected Array2DRowRealMatrix nordsieck;<a name="line.67"></a>
71 <FONT color="green">068</FONT> <a name="line.68"></a>
72 <FONT color="green">069</FONT> /** Starter integrator. */<a name="line.69"></a>
73 <FONT color="green">070</FONT> private FirstOrderIntegrator starter;<a name="line.70"></a>
74 <FONT color="green">071</FONT> <a name="line.71"></a>
75 <FONT color="green">072</FONT> /** Number of steps of the multistep method (excluding the one being computed). */<a name="line.72"></a>
76 <FONT color="green">073</FONT> private final int nSteps;<a name="line.73"></a>
77 <FONT color="green">074</FONT> <a name="line.74"></a>
78 <FONT color="green">075</FONT> /** Stepsize control exponent. */<a name="line.75"></a>
79 <FONT color="green">076</FONT> private double exp;<a name="line.76"></a>
80 <FONT color="green">077</FONT> <a name="line.77"></a>
81 <FONT color="green">078</FONT> /** Safety factor for stepsize control. */<a name="line.78"></a>
82 <FONT color="green">079</FONT> private double safety;<a name="line.79"></a>
83 <FONT color="green">080</FONT> <a name="line.80"></a>
84 <FONT color="green">081</FONT> /** Minimal reduction factor for stepsize control. */<a name="line.81"></a>
85 <FONT color="green">082</FONT> private double minReduction;<a name="line.82"></a>
86 <FONT color="green">083</FONT> <a name="line.83"></a>
87 <FONT color="green">084</FONT> /** Maximal growth factor for stepsize control. */<a name="line.84"></a>
88 <FONT color="green">085</FONT> private double maxGrowth;<a name="line.85"></a>
89 <FONT color="green">086</FONT> <a name="line.86"></a>
90 <FONT color="green">087</FONT> /**<a name="line.87"></a>
91 <FONT color="green">088</FONT> * Build a multistep integrator with the given stepsize bounds.<a name="line.88"></a>
92 <FONT color="green">089</FONT> * &lt;p&gt;The default starter integrator is set to the {@link<a name="line.89"></a>
93 <FONT color="green">090</FONT> * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with<a name="line.90"></a>
94 <FONT color="green">091</FONT> * some defaults settings.&lt;/p&gt;<a name="line.91"></a>
95 <FONT color="green">092</FONT> * &lt;p&gt;<a name="line.92"></a>
96 <FONT color="green">093</FONT> * The default max growth factor is set to a quite low value: 2&lt;sup&gt;1/order&lt;/sup&gt;.<a name="line.93"></a>
97 <FONT color="green">094</FONT> * &lt;/p&gt;<a name="line.94"></a>
98 <FONT color="green">095</FONT> * @param name name of the method<a name="line.95"></a>
99 <FONT color="green">096</FONT> * @param nSteps number of steps of the multistep method<a name="line.96"></a>
100 <FONT color="green">097</FONT> * (excluding the one being computed)<a name="line.97"></a>
101 <FONT color="green">098</FONT> * @param order order of the method<a name="line.98"></a>
102 <FONT color="green">099</FONT> * @param minStep minimal step (must be positive even for backward<a name="line.99"></a>
103 <FONT color="green">100</FONT> * integration), the last step can be smaller than this<a name="line.100"></a>
104 <FONT color="green">101</FONT> * @param maxStep maximal step (must be positive even for backward<a name="line.101"></a>
105 <FONT color="green">102</FONT> * integration)<a name="line.102"></a>
106 <FONT color="green">103</FONT> * @param scalAbsoluteTolerance allowed absolute error<a name="line.103"></a>
107 <FONT color="green">104</FONT> * @param scalRelativeTolerance allowed relative error<a name="line.104"></a>
108 <FONT color="green">105</FONT> */<a name="line.105"></a>
109 <FONT color="green">106</FONT> protected MultistepIntegrator(final String name, final int nSteps,<a name="line.106"></a>
110 <FONT color="green">107</FONT> final int order,<a name="line.107"></a>
111 <FONT color="green">108</FONT> final double minStep, final double maxStep,<a name="line.108"></a>
112 <FONT color="green">109</FONT> final double scalAbsoluteTolerance,<a name="line.109"></a>
113 <FONT color="green">110</FONT> final double scalRelativeTolerance) {<a name="line.110"></a>
114 <FONT color="green">111</FONT> <a name="line.111"></a>
115 <FONT color="green">112</FONT> super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);<a name="line.112"></a>
116 <FONT color="green">113</FONT> <a name="line.113"></a>
117 <FONT color="green">114</FONT> if (nSteps &lt;= 0) {<a name="line.114"></a>
118 <FONT color="green">115</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.115"></a>
119 <FONT color="green">116</FONT> "{0} method needs at least one previous point",<a name="line.116"></a>
120 <FONT color="green">117</FONT> name);<a name="line.117"></a>
121 <FONT color="green">118</FONT> }<a name="line.118"></a>
122 <FONT color="green">119</FONT> <a name="line.119"></a>
123 <FONT color="green">120</FONT> starter = new DormandPrince853Integrator(minStep, maxStep,<a name="line.120"></a>
124 <FONT color="green">121</FONT> scalAbsoluteTolerance,<a name="line.121"></a>
125 <FONT color="green">122</FONT> scalRelativeTolerance);<a name="line.122"></a>
126 <FONT color="green">123</FONT> this.nSteps = nSteps;<a name="line.123"></a>
127 <FONT color="green">124</FONT> <a name="line.124"></a>
128 <FONT color="green">125</FONT> exp = -1.0 / order;<a name="line.125"></a>
129 <FONT color="green">126</FONT> <a name="line.126"></a>
130 <FONT color="green">127</FONT> // set the default values of the algorithm control parameters<a name="line.127"></a>
131 <FONT color="green">128</FONT> setSafety(0.9);<a name="line.128"></a>
132 <FONT color="green">129</FONT> setMinReduction(0.2);<a name="line.129"></a>
133 <FONT color="green">130</FONT> setMaxGrowth(Math.pow(2.0, -exp));<a name="line.130"></a>
134 <FONT color="green">131</FONT> <a name="line.131"></a>
135 <FONT color="green">132</FONT> }<a name="line.132"></a>
136 <FONT color="green">133</FONT> <a name="line.133"></a>
137 <FONT color="green">134</FONT> /**<a name="line.134"></a>
138 <FONT color="green">135</FONT> * Build a multistep integrator with the given stepsize bounds.<a name="line.135"></a>
139 <FONT color="green">136</FONT> * &lt;p&gt;The default starter integrator is set to the {@link<a name="line.136"></a>
140 <FONT color="green">137</FONT> * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with<a name="line.137"></a>
141 <FONT color="green">138</FONT> * some defaults settings.&lt;/p&gt;<a name="line.138"></a>
142 <FONT color="green">139</FONT> * &lt;p&gt;<a name="line.139"></a>
143 <FONT color="green">140</FONT> * The default max growth factor is set to a quite low value: 2&lt;sup&gt;1/order&lt;/sup&gt;.<a name="line.140"></a>
144 <FONT color="green">141</FONT> * &lt;/p&gt;<a name="line.141"></a>
145 <FONT color="green">142</FONT> * @param name name of the method<a name="line.142"></a>
146 <FONT color="green">143</FONT> * @param nSteps number of steps of the multistep method<a name="line.143"></a>
147 <FONT color="green">144</FONT> * (excluding the one being computed)<a name="line.144"></a>
148 <FONT color="green">145</FONT> * @param order order of the method<a name="line.145"></a>
149 <FONT color="green">146</FONT> * @param minStep minimal step (must be positive even for backward<a name="line.146"></a>
150 <FONT color="green">147</FONT> * integration), the last step can be smaller than this<a name="line.147"></a>
151 <FONT color="green">148</FONT> * @param maxStep maximal step (must be positive even for backward<a name="line.148"></a>
152 <FONT color="green">149</FONT> * integration)<a name="line.149"></a>
153 <FONT color="green">150</FONT> * @param vecAbsoluteTolerance allowed absolute error<a name="line.150"></a>
154 <FONT color="green">151</FONT> * @param vecRelativeTolerance allowed relative error<a name="line.151"></a>
155 <FONT color="green">152</FONT> */<a name="line.152"></a>
156 <FONT color="green">153</FONT> protected MultistepIntegrator(final String name, final int nSteps,<a name="line.153"></a>
157 <FONT color="green">154</FONT> final int order,<a name="line.154"></a>
158 <FONT color="green">155</FONT> final double minStep, final double maxStep,<a name="line.155"></a>
159 <FONT color="green">156</FONT> final double[] vecAbsoluteTolerance,<a name="line.156"></a>
160 <FONT color="green">157</FONT> final double[] vecRelativeTolerance) {<a name="line.157"></a>
161 <FONT color="green">158</FONT> super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);<a name="line.158"></a>
162 <FONT color="green">159</FONT> starter = new DormandPrince853Integrator(minStep, maxStep,<a name="line.159"></a>
163 <FONT color="green">160</FONT> vecAbsoluteTolerance,<a name="line.160"></a>
164 <FONT color="green">161</FONT> vecRelativeTolerance);<a name="line.161"></a>
165 <FONT color="green">162</FONT> this.nSteps = nSteps;<a name="line.162"></a>
166 <FONT color="green">163</FONT> <a name="line.163"></a>
167 <FONT color="green">164</FONT> exp = -1.0 / order;<a name="line.164"></a>
168 <FONT color="green">165</FONT> <a name="line.165"></a>
169 <FONT color="green">166</FONT> // set the default values of the algorithm control parameters<a name="line.166"></a>
170 <FONT color="green">167</FONT> setSafety(0.9);<a name="line.167"></a>
171 <FONT color="green">168</FONT> setMinReduction(0.2);<a name="line.168"></a>
172 <FONT color="green">169</FONT> setMaxGrowth(Math.pow(2.0, -exp));<a name="line.169"></a>
173 <FONT color="green">170</FONT> <a name="line.170"></a>
174 <FONT color="green">171</FONT> }<a name="line.171"></a>
175 <FONT color="green">172</FONT> <a name="line.172"></a>
176 <FONT color="green">173</FONT> /**<a name="line.173"></a>
177 <FONT color="green">174</FONT> * Get the starter integrator.<a name="line.174"></a>
178 <FONT color="green">175</FONT> * @return starter integrator<a name="line.175"></a>
179 <FONT color="green">176</FONT> */<a name="line.176"></a>
180 <FONT color="green">177</FONT> public ODEIntegrator getStarterIntegrator() {<a name="line.177"></a>
181 <FONT color="green">178</FONT> return starter;<a name="line.178"></a>
182 <FONT color="green">179</FONT> }<a name="line.179"></a>
183 <FONT color="green">180</FONT> <a name="line.180"></a>
184 <FONT color="green">181</FONT> /**<a name="line.181"></a>
185 <FONT color="green">182</FONT> * Set the starter integrator.<a name="line.182"></a>
186 <FONT color="green">183</FONT> * &lt;p&gt;The various step and event handlers for this starter integrator<a name="line.183"></a>
187 <FONT color="green">184</FONT> * will be managed automatically by the multi-step integrator. Any<a name="line.184"></a>
188 <FONT color="green">185</FONT> * user configuration for these elements will be cleared before use.&lt;/p&gt;<a name="line.185"></a>
189 <FONT color="green">186</FONT> * @param starterIntegrator starter integrator<a name="line.186"></a>
190 <FONT color="green">187</FONT> */<a name="line.187"></a>
191 <FONT color="green">188</FONT> public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {<a name="line.188"></a>
192 <FONT color="green">189</FONT> this.starter = starterIntegrator;<a name="line.189"></a>
193 <FONT color="green">190</FONT> }<a name="line.190"></a>
194 <FONT color="green">191</FONT> <a name="line.191"></a>
195 <FONT color="green">192</FONT> /** Start the integration.<a name="line.192"></a>
196 <FONT color="green">193</FONT> * &lt;p&gt;This method computes one step using the underlying starter integrator,<a name="line.193"></a>
197 <FONT color="green">194</FONT> * and initializes the Nordsieck vector at step start. The starter integrator<a name="line.194"></a>
198 <FONT color="green">195</FONT> * purpose is only to establish initial conditions, it does not really change<a name="line.195"></a>
199 <FONT color="green">196</FONT> * time by itself. The top level multistep integrator remains in charge of<a name="line.196"></a>
200 <FONT color="green">197</FONT> * handling time propagation and events handling as it will starts its own<a name="line.197"></a>
201 <FONT color="green">198</FONT> * computation right from the beginning. In a sense, the starter integrator<a name="line.198"></a>
202 <FONT color="green">199</FONT> * can be seen as a dummy one and so it will never trigger any user event nor<a name="line.199"></a>
203 <FONT color="green">200</FONT> * call any user step handler.&lt;/p&gt;<a name="line.200"></a>
204 <FONT color="green">201</FONT> * @param t0 initial time<a name="line.201"></a>
205 <FONT color="green">202</FONT> * @param y0 initial value of the state vector at t0<a name="line.202"></a>
206 <FONT color="green">203</FONT> * @param t target time for the integration<a name="line.203"></a>
207 <FONT color="green">204</FONT> * (can be set to a value smaller than &lt;code&gt;t0&lt;/code&gt; for backward integration)<a name="line.204"></a>
208 <FONT color="green">205</FONT> * @throws IntegratorException if the integrator cannot perform integration<a name="line.205"></a>
209 <FONT color="green">206</FONT> * @throws DerivativeException this exception is propagated to the caller if<a name="line.206"></a>
210 <FONT color="green">207</FONT> * the underlying user function triggers one<a name="line.207"></a>
211 <FONT color="green">208</FONT> */<a name="line.208"></a>
212 <FONT color="green">209</FONT> protected void start(final double t0, final double[] y0, final double t)<a name="line.209"></a>
213 <FONT color="green">210</FONT> throws DerivativeException, IntegratorException {<a name="line.210"></a>
214 <FONT color="green">211</FONT> <a name="line.211"></a>
215 <FONT color="green">212</FONT> // make sure NO user event nor user step handler is triggered,<a name="line.212"></a>
216 <FONT color="green">213</FONT> // this is the task of the top level integrator, not the task<a name="line.213"></a>
217 <FONT color="green">214</FONT> // of the starter integrator<a name="line.214"></a>
218 <FONT color="green">215</FONT> starter.clearEventHandlers();<a name="line.215"></a>
219 <FONT color="green">216</FONT> starter.clearStepHandlers();<a name="line.216"></a>
220 <FONT color="green">217</FONT> <a name="line.217"></a>
221 <FONT color="green">218</FONT> // set up one specific step handler to extract initial Nordsieck vector<a name="line.218"></a>
222 <FONT color="green">219</FONT> starter.addStepHandler(new NordsieckInitializer(y0.length));<a name="line.219"></a>
223 <FONT color="green">220</FONT> <a name="line.220"></a>
224 <FONT color="green">221</FONT> // start integration, expecting a InitializationCompletedMarkerException<a name="line.221"></a>
225 <FONT color="green">222</FONT> try {<a name="line.222"></a>
226 <FONT color="green">223</FONT> starter.integrate(new CountingDifferentialEquations(y0.length),<a name="line.223"></a>
227 <FONT color="green">224</FONT> t0, y0, t, new double[y0.length]);<a name="line.224"></a>
228 <FONT color="green">225</FONT> } catch (DerivativeException de) {<a name="line.225"></a>
229 <FONT color="green">226</FONT> if (!(de instanceof InitializationCompletedMarkerException)) {<a name="line.226"></a>
230 <FONT color="green">227</FONT> // this is not the expected nominal interruption of the start integrator<a name="line.227"></a>
231 <FONT color="green">228</FONT> throw de;<a name="line.228"></a>
232 <FONT color="green">229</FONT> }<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> // remove the specific step handler<a name="line.232"></a>
236 <FONT color="green">233</FONT> starter.clearStepHandlers();<a name="line.233"></a>
237 <FONT color="green">234</FONT> <a name="line.234"></a>
238 <FONT color="green">235</FONT> }<a name="line.235"></a>
239 <FONT color="green">236</FONT> <a name="line.236"></a>
240 <FONT color="green">237</FONT> /** Initialize the high order scaled derivatives at step start.<a name="line.237"></a>
241 <FONT color="green">238</FONT> * @param first first scaled derivative at step start<a name="line.238"></a>
242 <FONT color="green">239</FONT> * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)<a name="line.239"></a>
243 <FONT color="green">240</FONT> * will be modified<a name="line.240"></a>
244 <FONT color="green">241</FONT> * @return high order scaled derivatives at step start<a name="line.241"></a>
245 <FONT color="green">242</FONT> */<a name="line.242"></a>
246 <FONT color="green">243</FONT> protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,<a name="line.243"></a>
247 <FONT color="green">244</FONT> final double[][] multistep);<a name="line.244"></a>
248 <FONT color="green">245</FONT> <a name="line.245"></a>
249 <FONT color="green">246</FONT> /** Get the minimal reduction factor for stepsize control.<a name="line.246"></a>
250 <FONT color="green">247</FONT> * @return minimal reduction factor<a name="line.247"></a>
251 <FONT color="green">248</FONT> */<a name="line.248"></a>
252 <FONT color="green">249</FONT> public double getMinReduction() {<a name="line.249"></a>
253 <FONT color="green">250</FONT> return minReduction;<a name="line.250"></a>
254 <FONT color="green">251</FONT> }<a name="line.251"></a>
255 <FONT color="green">252</FONT> <a name="line.252"></a>
256 <FONT color="green">253</FONT> /** Set the minimal reduction factor for stepsize control.<a name="line.253"></a>
257 <FONT color="green">254</FONT> * @param minReduction minimal reduction factor<a name="line.254"></a>
258 <FONT color="green">255</FONT> */<a name="line.255"></a>
259 <FONT color="green">256</FONT> public void setMinReduction(final double minReduction) {<a name="line.256"></a>
260 <FONT color="green">257</FONT> this.minReduction = minReduction;<a name="line.257"></a>
261 <FONT color="green">258</FONT> }<a name="line.258"></a>
262 <FONT color="green">259</FONT> <a name="line.259"></a>
263 <FONT color="green">260</FONT> /** Get the maximal growth factor for stepsize control.<a name="line.260"></a>
264 <FONT color="green">261</FONT> * @return maximal growth factor<a name="line.261"></a>
265 <FONT color="green">262</FONT> */<a name="line.262"></a>
266 <FONT color="green">263</FONT> public double getMaxGrowth() {<a name="line.263"></a>
267 <FONT color="green">264</FONT> return maxGrowth;<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> /** Set the maximal growth factor for stepsize control.<a name="line.267"></a>
271 <FONT color="green">268</FONT> * @param maxGrowth maximal growth factor<a name="line.268"></a>
272 <FONT color="green">269</FONT> */<a name="line.269"></a>
273 <FONT color="green">270</FONT> public void setMaxGrowth(final double maxGrowth) {<a name="line.270"></a>
274 <FONT color="green">271</FONT> this.maxGrowth = maxGrowth;<a name="line.271"></a>
275 <FONT color="green">272</FONT> }<a name="line.272"></a>
276 <FONT color="green">273</FONT> <a name="line.273"></a>
277 <FONT color="green">274</FONT> /** Get the safety factor for stepsize control.<a name="line.274"></a>
278 <FONT color="green">275</FONT> * @return safety factor<a name="line.275"></a>
279 <FONT color="green">276</FONT> */<a name="line.276"></a>
280 <FONT color="green">277</FONT> public double getSafety() {<a name="line.277"></a>
281 <FONT color="green">278</FONT> return safety;<a name="line.278"></a>
282 <FONT color="green">279</FONT> }<a name="line.279"></a>
283 <FONT color="green">280</FONT> <a name="line.280"></a>
284 <FONT color="green">281</FONT> /** Set the safety factor for stepsize control.<a name="line.281"></a>
285 <FONT color="green">282</FONT> * @param safety safety factor<a name="line.282"></a>
286 <FONT color="green">283</FONT> */<a name="line.283"></a>
287 <FONT color="green">284</FONT> public void setSafety(final double safety) {<a name="line.284"></a>
288 <FONT color="green">285</FONT> this.safety = safety;<a name="line.285"></a>
289 <FONT color="green">286</FONT> }<a name="line.286"></a>
290 <FONT color="green">287</FONT> <a name="line.287"></a>
291 <FONT color="green">288</FONT> /** Compute step grow/shrink factor according to normalized error.<a name="line.288"></a>
292 <FONT color="green">289</FONT> * @param error normalized error of the current step<a name="line.289"></a>
293 <FONT color="green">290</FONT> * @return grow/shrink factor for next step<a name="line.290"></a>
294 <FONT color="green">291</FONT> */<a name="line.291"></a>
295 <FONT color="green">292</FONT> protected double computeStepGrowShrinkFactor(final double error) {<a name="line.292"></a>
296 <FONT color="green">293</FONT> return Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp)));<a name="line.293"></a>
297 <FONT color="green">294</FONT> }<a name="line.294"></a>
298 <FONT color="green">295</FONT> <a name="line.295"></a>
299 <FONT color="green">296</FONT> /** Transformer used to convert the first step to Nordsieck representation. */<a name="line.296"></a>
300 <FONT color="green">297</FONT> public static interface NordsieckTransformer {<a name="line.297"></a>
301 <FONT color="green">298</FONT> /** Initialize the high order scaled derivatives at step start.<a name="line.298"></a>
302 <FONT color="green">299</FONT> * @param first first scaled derivative at step start<a name="line.299"></a>
303 <FONT color="green">300</FONT> * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)<a name="line.300"></a>
304 <FONT color="green">301</FONT> * will be modified<a name="line.301"></a>
305 <FONT color="green">302</FONT> * @return high order derivatives at step start<a name="line.302"></a>
306 <FONT color="green">303</FONT> */<a name="line.303"></a>
307 <FONT color="green">304</FONT> RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);<a name="line.304"></a>
308 <FONT color="green">305</FONT> }<a name="line.305"></a>
309 <FONT color="green">306</FONT> <a name="line.306"></a>
310 <FONT color="green">307</FONT> /** Specialized step handler storing the first step. */<a name="line.307"></a>
311 <FONT color="green">308</FONT> private class NordsieckInitializer implements StepHandler {<a name="line.308"></a>
312 <FONT color="green">309</FONT> <a name="line.309"></a>
313 <FONT color="green">310</FONT> /** Problem dimension. */<a name="line.310"></a>
314 <FONT color="green">311</FONT> private final int n;<a name="line.311"></a>
315 <FONT color="green">312</FONT> <a name="line.312"></a>
316 <FONT color="green">313</FONT> /** Simple constructor.<a name="line.313"></a>
317 <FONT color="green">314</FONT> * @param n problem dimension<a name="line.314"></a>
318 <FONT color="green">315</FONT> */<a name="line.315"></a>
319 <FONT color="green">316</FONT> public NordsieckInitializer(final int n) {<a name="line.316"></a>
320 <FONT color="green">317</FONT> this.n = n;<a name="line.317"></a>
321 <FONT color="green">318</FONT> }<a name="line.318"></a>
322 <FONT color="green">319</FONT> <a name="line.319"></a>
323 <FONT color="green">320</FONT> /** {@inheritDoc} */<a name="line.320"></a>
324 <FONT color="green">321</FONT> public void handleStep(StepInterpolator interpolator, boolean isLast)<a name="line.321"></a>
325 <FONT color="green">322</FONT> throws DerivativeException {<a name="line.322"></a>
326 <FONT color="green">323</FONT> <a name="line.323"></a>
327 <FONT color="green">324</FONT> final double prev = interpolator.getPreviousTime();<a name="line.324"></a>
328 <FONT color="green">325</FONT> final double curr = interpolator.getCurrentTime();<a name="line.325"></a>
329 <FONT color="green">326</FONT> stepStart = prev;<a name="line.326"></a>
330 <FONT color="green">327</FONT> stepSize = (curr - prev) / (nSteps + 1);<a name="line.327"></a>
331 <FONT color="green">328</FONT> <a name="line.328"></a>
332 <FONT color="green">329</FONT> // compute the first scaled derivative<a name="line.329"></a>
333 <FONT color="green">330</FONT> interpolator.setInterpolatedTime(prev);<a name="line.330"></a>
334 <FONT color="green">331</FONT> scaled = interpolator.getInterpolatedDerivatives().clone();<a name="line.331"></a>
335 <FONT color="green">332</FONT> for (int j = 0; j &lt; n; ++j) {<a name="line.332"></a>
336 <FONT color="green">333</FONT> scaled[j] *= stepSize;<a name="line.333"></a>
337 <FONT color="green">334</FONT> }<a name="line.334"></a>
338 <FONT color="green">335</FONT> <a name="line.335"></a>
339 <FONT color="green">336</FONT> // compute the high order scaled derivatives<a name="line.336"></a>
340 <FONT color="green">337</FONT> final double[][] multistep = new double[nSteps][];<a name="line.337"></a>
341 <FONT color="green">338</FONT> for (int i = 1; i &lt;= nSteps; ++i) {<a name="line.338"></a>
342 <FONT color="green">339</FONT> interpolator.setInterpolatedTime(prev + stepSize * i);<a name="line.339"></a>
343 <FONT color="green">340</FONT> final double[] msI = interpolator.getInterpolatedDerivatives().clone();<a name="line.340"></a>
344 <FONT color="green">341</FONT> for (int j = 0; j &lt; n; ++j) {<a name="line.341"></a>
345 <FONT color="green">342</FONT> msI[j] *= stepSize;<a name="line.342"></a>
346 <FONT color="green">343</FONT> }<a name="line.343"></a>
347 <FONT color="green">344</FONT> multistep[i - 1] = msI;<a name="line.344"></a>
348 <FONT color="green">345</FONT> }<a name="line.345"></a>
349 <FONT color="green">346</FONT> nordsieck = initializeHighOrderDerivatives(scaled, multistep);<a name="line.346"></a>
350 <FONT color="green">347</FONT> <a name="line.347"></a>
351 <FONT color="green">348</FONT> // stop the integrator after the first step has been handled<a name="line.348"></a>
352 <FONT color="green">349</FONT> throw new InitializationCompletedMarkerException();<a name="line.349"></a>
353 <FONT color="green">350</FONT> <a name="line.350"></a>
354 <FONT color="green">351</FONT> }<a name="line.351"></a>
355 <FONT color="green">352</FONT> <a name="line.352"></a>
356 <FONT color="green">353</FONT> /** {@inheritDoc} */<a name="line.353"></a>
357 <FONT color="green">354</FONT> public boolean requiresDenseOutput() {<a name="line.354"></a>
358 <FONT color="green">355</FONT> return true;<a name="line.355"></a>
359 <FONT color="green">356</FONT> }<a name="line.356"></a>
360 <FONT color="green">357</FONT> <a name="line.357"></a>
361 <FONT color="green">358</FONT> /** {@inheritDoc} */<a name="line.358"></a>
362 <FONT color="green">359</FONT> public void reset() {<a name="line.359"></a>
363 <FONT color="green">360</FONT> // nothing to do<a name="line.360"></a>
364 <FONT color="green">361</FONT> }<a name="line.361"></a>
365 <FONT color="green">362</FONT> <a name="line.362"></a>
366 <FONT color="green">363</FONT> }<a name="line.363"></a>
367 <FONT color="green">364</FONT> <a name="line.364"></a>
368 <FONT color="green">365</FONT> /** Marker exception used ONLY to stop the starter integrator after first step. */<a name="line.365"></a>
369 <FONT color="green">366</FONT> private static class InitializationCompletedMarkerException<a name="line.366"></a>
370 <FONT color="green">367</FONT> extends DerivativeException {<a name="line.367"></a>
371 <FONT color="green">368</FONT> <a name="line.368"></a>
372 <FONT color="green">369</FONT> /** Serializable version identifier. */<a name="line.369"></a>
373 <FONT color="green">370</FONT> private static final long serialVersionUID = -4105805787353488365L;<a name="line.370"></a>
374 <FONT color="green">371</FONT> <a name="line.371"></a>
375 <FONT color="green">372</FONT> /** Simple constructor. */<a name="line.372"></a>
376 <FONT color="green">373</FONT> public InitializationCompletedMarkerException() {<a name="line.373"></a>
377 <FONT color="green">374</FONT> super((Throwable) null);<a name="line.374"></a>
378 <FONT color="green">375</FONT> }<a name="line.375"></a>
379 <FONT color="green">376</FONT> <a name="line.376"></a>
380 <FONT color="green">377</FONT> }<a name="line.377"></a>
381 <FONT color="green">378</FONT> <a name="line.378"></a>
382 <FONT color="green">379</FONT> /** Wrapper for differential equations, ensuring start evaluations are counted. */<a name="line.379"></a>
383 <FONT color="green">380</FONT> private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {<a name="line.380"></a>
384 <FONT color="green">381</FONT> <a name="line.381"></a>
385 <FONT color="green">382</FONT> /** Dimension of the problem. */<a name="line.382"></a>
386 <FONT color="green">383</FONT> private final int dimension;<a name="line.383"></a>
387 <FONT color="green">384</FONT> <a name="line.384"></a>
388 <FONT color="green">385</FONT> /** Simple constructor.<a name="line.385"></a>
389 <FONT color="green">386</FONT> * @param dimension dimension of the problem<a name="line.386"></a>
390 <FONT color="green">387</FONT> */<a name="line.387"></a>
391 <FONT color="green">388</FONT> public CountingDifferentialEquations(final int dimension) {<a name="line.388"></a>
392 <FONT color="green">389</FONT> this.dimension = dimension;<a name="line.389"></a>
393 <FONT color="green">390</FONT> }<a name="line.390"></a>
394 <FONT color="green">391</FONT> <a name="line.391"></a>
395 <FONT color="green">392</FONT> /** {@inheritDoc} */<a name="line.392"></a>
396 <FONT color="green">393</FONT> public void computeDerivatives(double t, double[] y, double[] dot)<a name="line.393"></a>
397 <FONT color="green">394</FONT> throws DerivativeException {<a name="line.394"></a>
398 <FONT color="green">395</FONT> MultistepIntegrator.this.computeDerivatives(t, y, dot);<a name="line.395"></a>
399 <FONT color="green">396</FONT> }<a name="line.396"></a>
400 <FONT color="green">397</FONT> <a name="line.397"></a>
401 <FONT color="green">398</FONT> /** {@inheritDoc} */<a name="line.398"></a>
402 <FONT color="green">399</FONT> public int getDimension() {<a name="line.399"></a>
403 <FONT color="green">400</FONT> return dimension;<a name="line.400"></a>
404 <FONT color="green">401</FONT> }<a name="line.401"></a>
405 <FONT color="green">402</FONT> }<a name="line.402"></a>
406 <FONT color="green">403</FONT> <a name="line.403"></a>
407 <FONT color="green">404</FONT> }<a name="line.404"></a>
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468 </PRE>
469 </BODY>
470 </HTML>