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

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.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.linear.Array2DRowRealMatrix;<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> import org.apache.commons.math.ode.events.CombinedEventsManager;<a name="line.24"></a>
28 <FONT color="green">025</FONT> import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;<a name="line.25"></a>
29 <FONT color="green">026</FONT> import org.apache.commons.math.ode.sampling.StepHandler;<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> /**<a name="line.29"></a>
33 <FONT color="green">030</FONT> * This class implements explicit Adams-Bashforth integrators for Ordinary<a name="line.30"></a>
34 <FONT color="green">031</FONT> * Differential Equations.<a name="line.31"></a>
35 <FONT color="green">032</FONT> *<a name="line.32"></a>
36 <FONT color="green">033</FONT> * &lt;p&gt;Adams-Bashforth methods (in fact due to Adams alone) are explicit<a name="line.33"></a>
37 <FONT color="green">034</FONT> * multistep ODE solvers. This implementation is a variation of the classical<a name="line.34"></a>
38 <FONT color="green">035</FONT> * one: it uses adaptive stepsize to implement error control, whereas<a name="line.35"></a>
39 <FONT color="green">036</FONT> * classical implementations are fixed step size. The value of state vector<a name="line.36"></a>
40 <FONT color="green">037</FONT> * at step n+1 is a simple combination of the value at step n and of the<a name="line.37"></a>
41 <FONT color="green">038</FONT> * derivatives at steps n, n-1, n-2 ... Depending on the number k of previous<a name="line.38"></a>
42 <FONT color="green">039</FONT> * steps one wants to use for computing the next value, different formulas<a name="line.39"></a>
43 <FONT color="green">040</FONT> * are available:&lt;/p&gt;<a name="line.40"></a>
44 <FONT color="green">041</FONT> * &lt;ul&gt;<a name="line.41"></a>
45 <FONT color="green">042</FONT> * &lt;li&gt;k = 1: y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + h y'&lt;sub&gt;n&lt;/sub&gt;&lt;/li&gt;<a name="line.42"></a>
46 <FONT color="green">043</FONT> * &lt;li&gt;k = 2: y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + h (3y'&lt;sub&gt;n&lt;/sub&gt;-y'&lt;sub&gt;n-1&lt;/sub&gt;)/2&lt;/li&gt;<a name="line.43"></a>
47 <FONT color="green">044</FONT> * &lt;li&gt;k = 3: y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + h (23y'&lt;sub&gt;n&lt;/sub&gt;-16y'&lt;sub&gt;n-1&lt;/sub&gt;+5y'&lt;sub&gt;n-2&lt;/sub&gt;)/12&lt;/li&gt;<a name="line.44"></a>
48 <FONT color="green">045</FONT> * &lt;li&gt;k = 4: y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + h (55y'&lt;sub&gt;n&lt;/sub&gt;-59y'&lt;sub&gt;n-1&lt;/sub&gt;+37y'&lt;sub&gt;n-2&lt;/sub&gt;-9y'&lt;sub&gt;n-3&lt;/sub&gt;)/24&lt;/li&gt;<a name="line.45"></a>
49 <FONT color="green">046</FONT> * &lt;li&gt;...&lt;/li&gt;<a name="line.46"></a>
50 <FONT color="green">047</FONT> * &lt;/ul&gt;<a name="line.47"></a>
51 <FONT color="green">048</FONT> *<a name="line.48"></a>
52 <FONT color="green">049</FONT> * &lt;p&gt;A k-steps Adams-Bashforth method is of order k.&lt;/p&gt;<a name="line.49"></a>
53 <FONT color="green">050</FONT> *<a name="line.50"></a>
54 <FONT color="green">051</FONT> * &lt;h3&gt;Implementation details&lt;/h3&gt;<a name="line.51"></a>
55 <FONT color="green">052</FONT> *<a name="line.52"></a>
56 <FONT color="green">053</FONT> * &lt;p&gt;We define scaled derivatives s&lt;sub&gt;i&lt;/sub&gt;(n) at step n as:<a name="line.53"></a>
57 <FONT color="green">054</FONT> * &lt;pre&gt;<a name="line.54"></a>
58 <FONT color="green">055</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.55"></a>
59 <FONT color="green">056</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.56"></a>
60 <FONT color="green">057</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.57"></a>
61 <FONT color="green">058</FONT> * ...<a name="line.58"></a>
62 <FONT color="green">059</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.59"></a>
63 <FONT color="green">060</FONT> * &lt;/pre&gt;&lt;/p&gt;<a name="line.60"></a>
64 <FONT color="green">061</FONT> *<a name="line.61"></a>
65 <FONT color="green">062</FONT> * &lt;p&gt;The definitions above use the classical representation with several previous first<a name="line.62"></a>
66 <FONT color="green">063</FONT> * derivatives. Lets define<a name="line.63"></a>
67 <FONT color="green">064</FONT> * &lt;pre&gt;<a name="line.64"></a>
68 <FONT color="green">065</FONT> * q&lt;sub&gt;n&lt;/sub&gt; = [ s&lt;sub&gt;1&lt;/sub&gt;(n-1) s&lt;sub&gt;1&lt;/sub&gt;(n-2) ... s&lt;sub&gt;1&lt;/sub&gt;(n-(k-1)) ]&lt;sup&gt;T&lt;/sup&gt;<a name="line.65"></a>
69 <FONT color="green">066</FONT> * &lt;/pre&gt;<a name="line.66"></a>
70 <FONT color="green">067</FONT> * (we omit the k index in the notation for clarity). With these definitions,<a name="line.67"></a>
71 <FONT color="green">068</FONT> * Adams-Bashforth methods can be written:<a name="line.68"></a>
72 <FONT color="green">069</FONT> * &lt;ul&gt;<a name="line.69"></a>
73 <FONT color="green">070</FONT> * &lt;li&gt;k = 1: y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + s&lt;sub&gt;1&lt;/sub&gt;(n)&lt;/li&gt;<a name="line.70"></a>
74 <FONT color="green">071</FONT> * &lt;li&gt;k = 2: y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + 3/2 s&lt;sub&gt;1&lt;/sub&gt;(n) + [ -1/2 ] q&lt;sub&gt;n&lt;/sub&gt;&lt;/li&gt;<a name="line.71"></a>
75 <FONT color="green">072</FONT> * &lt;li&gt;k = 3: y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + 23/12 s&lt;sub&gt;1&lt;/sub&gt;(n) + [ -16/12 5/12 ] q&lt;sub&gt;n&lt;/sub&gt;&lt;/li&gt;<a name="line.72"></a>
76 <FONT color="green">073</FONT> * &lt;li&gt;k = 4: y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + 55/24 s&lt;sub&gt;1&lt;/sub&gt;(n) + [ -59/24 37/24 -9/24 ] q&lt;sub&gt;n&lt;/sub&gt;&lt;/li&gt;<a name="line.73"></a>
77 <FONT color="green">074</FONT> * &lt;li&gt;...&lt;/li&gt;<a name="line.74"></a>
78 <FONT color="green">075</FONT> * &lt;/ul&gt;&lt;/p&gt;<a name="line.75"></a>
79 <FONT color="green">076</FONT> *<a name="line.76"></a>
80 <FONT color="green">077</FONT> * &lt;p&gt;Instead of using the classical representation with first derivatives only (y&lt;sub&gt;n&lt;/sub&gt;,<a name="line.77"></a>
81 <FONT color="green">078</FONT> * s&lt;sub&gt;1&lt;/sub&gt;(n) and q&lt;sub&gt;n&lt;/sub&gt;), our implementation uses the Nordsieck vector with<a name="line.78"></a>
82 <FONT color="green">079</FONT> * higher degrees scaled derivatives all taken at the same step (y&lt;sub&gt;n&lt;/sub&gt;, s&lt;sub&gt;1&lt;/sub&gt;(n)<a name="line.79"></a>
83 <FONT color="green">080</FONT> * and r&lt;sub&gt;n&lt;/sub&gt;) where r&lt;sub&gt;n&lt;/sub&gt; is defined as:<a name="line.80"></a>
84 <FONT color="green">081</FONT> * &lt;pre&gt;<a name="line.81"></a>
85 <FONT color="green">082</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.82"></a>
86 <FONT color="green">083</FONT> * &lt;/pre&gt;<a name="line.83"></a>
87 <FONT color="green">084</FONT> * (here again we omit the k index in the notation for clarity)<a name="line.84"></a>
88 <FONT color="green">085</FONT> * &lt;/p&gt;<a name="line.85"></a>
89 <FONT color="green">086</FONT> *<a name="line.86"></a>
90 <FONT color="green">087</FONT> * &lt;p&gt;Taylor series formulas show that for any index offset i, s&lt;sub&gt;1&lt;/sub&gt;(n-i) can be<a name="line.87"></a>
91 <FONT color="green">088</FONT> * computed from s&lt;sub&gt;1&lt;/sub&gt;(n), s&lt;sub&gt;2&lt;/sub&gt;(n) ... s&lt;sub&gt;k&lt;/sub&gt;(n), the formula being exact<a name="line.88"></a>
92 <FONT color="green">089</FONT> * for degree k polynomials.<a name="line.89"></a>
93 <FONT color="green">090</FONT> * &lt;pre&gt;<a name="line.90"></a>
94 <FONT color="green">091</FONT> * s&lt;sub&gt;1&lt;/sub&gt;(n-i) = s&lt;sub&gt;1&lt;/sub&gt;(n) + &amp;sum;&lt;sub&gt;j&lt;/sub&gt; j (-i)&lt;sup&gt;j-1&lt;/sup&gt; s&lt;sub&gt;j&lt;/sub&gt;(n)<a name="line.91"></a>
95 <FONT color="green">092</FONT> * &lt;/pre&gt;<a name="line.92"></a>
96 <FONT color="green">093</FONT> * The previous formula can be used with several values for i to compute the transform between<a name="line.93"></a>
97 <FONT color="green">094</FONT> * classical representation and Nordsieck vector. The transform between r&lt;sub&gt;n&lt;/sub&gt;<a name="line.94"></a>
98 <FONT color="green">095</FONT> * and q&lt;sub&gt;n&lt;/sub&gt; resulting from the Taylor series formulas above is:<a name="line.95"></a>
99 <FONT color="green">096</FONT> * &lt;pre&gt;<a name="line.96"></a>
100 <FONT color="green">097</FONT> * q&lt;sub&gt;n&lt;/sub&gt; = s&lt;sub&gt;1&lt;/sub&gt;(n) u + P r&lt;sub&gt;n&lt;/sub&gt;<a name="line.97"></a>
101 <FONT color="green">098</FONT> * &lt;/pre&gt;<a name="line.98"></a>
102 <FONT color="green">099</FONT> * where u is the [ 1 1 ... 1 ]&lt;sup&gt;T&lt;/sup&gt; vector and P is the (k-1)&amp;times;(k-1) matrix built<a name="line.99"></a>
103 <FONT color="green">100</FONT> * with the j (-i)&lt;sup&gt;j-1&lt;/sup&gt; terms:<a name="line.100"></a>
104 <FONT color="green">101</FONT> * &lt;pre&gt;<a name="line.101"></a>
105 <FONT color="green">102</FONT> * [ -2 3 -4 5 ... ]<a name="line.102"></a>
106 <FONT color="green">103</FONT> * [ -4 12 -32 80 ... ]<a name="line.103"></a>
107 <FONT color="green">104</FONT> * P = [ -6 27 -108 405 ... ]<a name="line.104"></a>
108 <FONT color="green">105</FONT> * [ -8 48 -256 1280 ... ]<a name="line.105"></a>
109 <FONT color="green">106</FONT> * [ ... ]<a name="line.106"></a>
110 <FONT color="green">107</FONT> * &lt;/pre&gt;&lt;/p&gt;<a name="line.107"></a>
111 <FONT color="green">108</FONT> *<a name="line.108"></a>
112 <FONT color="green">109</FONT> * &lt;p&gt;Using the Nordsieck vector has several advantages:<a name="line.109"></a>
113 <FONT color="green">110</FONT> * &lt;ul&gt;<a name="line.110"></a>
114 <FONT color="green">111</FONT> * &lt;li&gt;it greatly simplifies step interpolation as the interpolator mainly applies<a name="line.111"></a>
115 <FONT color="green">112</FONT> * Taylor series formulas,&lt;/li&gt;<a name="line.112"></a>
116 <FONT color="green">113</FONT> * &lt;li&gt;it simplifies step changes that occur when discrete events that truncate<a name="line.113"></a>
117 <FONT color="green">114</FONT> * the step are triggered,&lt;/li&gt;<a name="line.114"></a>
118 <FONT color="green">115</FONT> * &lt;li&gt;it allows to extend the methods in order to support adaptive stepsize.&lt;/li&gt;<a name="line.115"></a>
119 <FONT color="green">116</FONT> * &lt;/ul&gt;&lt;/p&gt;<a name="line.116"></a>
120 <FONT color="green">117</FONT> *<a name="line.117"></a>
121 <FONT color="green">118</FONT> * &lt;p&gt;The Nordsieck vector at step n+1 is computed from the Nordsieck vector at step n as follows:<a name="line.118"></a>
122 <FONT color="green">119</FONT> * &lt;ul&gt;<a name="line.119"></a>
123 <FONT color="green">120</FONT> * &lt;li&gt;y&lt;sub&gt;n+1&lt;/sub&gt; = y&lt;sub&gt;n&lt;/sub&gt; + s&lt;sub&gt;1&lt;/sub&gt;(n) + u&lt;sup&gt;T&lt;/sup&gt; r&lt;sub&gt;n&lt;/sub&gt;&lt;/li&gt;<a name="line.120"></a>
124 <FONT color="green">121</FONT> * &lt;li&gt;s&lt;sub&gt;1&lt;/sub&gt;(n+1) = h f(t&lt;sub&gt;n+1&lt;/sub&gt;, y&lt;sub&gt;n+1&lt;/sub&gt;)&lt;/li&gt;<a name="line.121"></a>
125 <FONT color="green">122</FONT> * &lt;li&gt;r&lt;sub&gt;n+1&lt;/sub&gt; = (s&lt;sub&gt;1&lt;/sub&gt;(n) - s&lt;sub&gt;1&lt;/sub&gt;(n+1)) P&lt;sup&gt;-1&lt;/sup&gt; u + P&lt;sup&gt;-1&lt;/sup&gt; A P r&lt;sub&gt;n&lt;/sub&gt;&lt;/li&gt;<a name="line.122"></a>
126 <FONT color="green">123</FONT> * &lt;/ul&gt;<a name="line.123"></a>
127 <FONT color="green">124</FONT> * where A is a rows shifting matrix (the lower left part is an identity matrix):<a name="line.124"></a>
128 <FONT color="green">125</FONT> * &lt;pre&gt;<a name="line.125"></a>
129 <FONT color="green">126</FONT> * [ 0 0 ... 0 0 | 0 ]<a name="line.126"></a>
130 <FONT color="green">127</FONT> * [ ---------------+---]<a name="line.127"></a>
131 <FONT color="green">128</FONT> * [ 1 0 ... 0 0 | 0 ]<a name="line.128"></a>
132 <FONT color="green">129</FONT> * A = [ 0 1 ... 0 0 | 0 ]<a name="line.129"></a>
133 <FONT color="green">130</FONT> * [ ... | 0 ]<a name="line.130"></a>
134 <FONT color="green">131</FONT> * [ 0 0 ... 1 0 | 0 ]<a name="line.131"></a>
135 <FONT color="green">132</FONT> * [ 0 0 ... 0 1 | 0 ]<a name="line.132"></a>
136 <FONT color="green">133</FONT> * &lt;/pre&gt;&lt;/p&gt;<a name="line.133"></a>
137 <FONT color="green">134</FONT> *<a name="line.134"></a>
138 <FONT color="green">135</FONT> * &lt;p&gt;The P&lt;sup&gt;-1&lt;/sup&gt;u vector and the P&lt;sup&gt;-1&lt;/sup&gt; A P matrix do not depend on the state,<a name="line.135"></a>
139 <FONT color="green">136</FONT> * they only depend on k and therefore are precomputed once for all.&lt;/p&gt;<a name="line.136"></a>
140 <FONT color="green">137</FONT> *<a name="line.137"></a>
141 <FONT color="green">138</FONT> * @version $Revision: 927202 $ $Date: 2010-03-24 18:11:51 -0400 (Wed, 24 Mar 2010) $<a name="line.138"></a>
142 <FONT color="green">139</FONT> * @since 2.0<a name="line.139"></a>
143 <FONT color="green">140</FONT> */<a name="line.140"></a>
144 <FONT color="green">141</FONT> public class AdamsBashforthIntegrator extends AdamsIntegrator {<a name="line.141"></a>
145 <FONT color="green">142</FONT> <a name="line.142"></a>
146 <FONT color="green">143</FONT> /**<a name="line.143"></a>
147 <FONT color="green">144</FONT> * Build an Adams-Bashforth integrator with the given order and step control parameters.<a name="line.144"></a>
148 <FONT color="green">145</FONT> * @param nSteps number of steps of the method excluding the one being computed<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 scalAbsoluteTolerance allowed absolute error<a name="line.150"></a>
154 <FONT color="green">151</FONT> * @param scalRelativeTolerance allowed relative error<a name="line.151"></a>
155 <FONT color="green">152</FONT> * @exception IllegalArgumentException if order is 1 or less<a name="line.152"></a>
156 <FONT color="green">153</FONT> */<a name="line.153"></a>
157 <FONT color="green">154</FONT> public AdamsBashforthIntegrator(final int nSteps,<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 scalAbsoluteTolerance,<a name="line.156"></a>
160 <FONT color="green">157</FONT> final double scalRelativeTolerance)<a name="line.157"></a>
161 <FONT color="green">158</FONT> throws IllegalArgumentException {<a name="line.158"></a>
162 <FONT color="green">159</FONT> super("Adams-Bashforth", nSteps, nSteps, minStep, maxStep,<a name="line.159"></a>
163 <FONT color="green">160</FONT> scalAbsoluteTolerance, scalRelativeTolerance);<a name="line.160"></a>
164 <FONT color="green">161</FONT> }<a name="line.161"></a>
165 <FONT color="green">162</FONT> <a name="line.162"></a>
166 <FONT color="green">163</FONT> /**<a name="line.163"></a>
167 <FONT color="green">164</FONT> * Build an Adams-Bashforth integrator with the given order and step control parameters.<a name="line.164"></a>
168 <FONT color="green">165</FONT> * @param nSteps number of steps of the method excluding the one being computed<a name="line.165"></a>
169 <FONT color="green">166</FONT> * @param minStep minimal step (must be positive even for backward<a name="line.166"></a>
170 <FONT color="green">167</FONT> * integration), the last step can be smaller than this<a name="line.167"></a>
171 <FONT color="green">168</FONT> * @param maxStep maximal step (must be positive even for backward<a name="line.168"></a>
172 <FONT color="green">169</FONT> * integration)<a name="line.169"></a>
173 <FONT color="green">170</FONT> * @param vecAbsoluteTolerance allowed absolute error<a name="line.170"></a>
174 <FONT color="green">171</FONT> * @param vecRelativeTolerance allowed relative error<a name="line.171"></a>
175 <FONT color="green">172</FONT> * @exception IllegalArgumentException if order is 1 or less<a name="line.172"></a>
176 <FONT color="green">173</FONT> */<a name="line.173"></a>
177 <FONT color="green">174</FONT> public AdamsBashforthIntegrator(final int nSteps,<a name="line.174"></a>
178 <FONT color="green">175</FONT> final double minStep, final double maxStep,<a name="line.175"></a>
179 <FONT color="green">176</FONT> final double[] vecAbsoluteTolerance,<a name="line.176"></a>
180 <FONT color="green">177</FONT> final double[] vecRelativeTolerance)<a name="line.177"></a>
181 <FONT color="green">178</FONT> throws IllegalArgumentException {<a name="line.178"></a>
182 <FONT color="green">179</FONT> super("Adams-Bashforth", nSteps, nSteps, minStep, maxStep,<a name="line.179"></a>
183 <FONT color="green">180</FONT> vecAbsoluteTolerance, vecRelativeTolerance);<a name="line.180"></a>
184 <FONT color="green">181</FONT> }<a name="line.181"></a>
185 <FONT color="green">182</FONT> <a name="line.182"></a>
186 <FONT color="green">183</FONT> /** {@inheritDoc} */<a name="line.183"></a>
187 <FONT color="green">184</FONT> @Override<a name="line.184"></a>
188 <FONT color="green">185</FONT> public double integrate(final FirstOrderDifferentialEquations equations,<a name="line.185"></a>
189 <FONT color="green">186</FONT> final double t0, final double[] y0,<a name="line.186"></a>
190 <FONT color="green">187</FONT> final double t, final double[] y)<a name="line.187"></a>
191 <FONT color="green">188</FONT> throws DerivativeException, IntegratorException {<a name="line.188"></a>
192 <FONT color="green">189</FONT> <a name="line.189"></a>
193 <FONT color="green">190</FONT> final int n = y0.length;<a name="line.190"></a>
194 <FONT color="green">191</FONT> sanityChecks(equations, t0, y0, t, y);<a name="line.191"></a>
195 <FONT color="green">192</FONT> setEquations(equations);<a name="line.192"></a>
196 <FONT color="green">193</FONT> resetEvaluations();<a name="line.193"></a>
197 <FONT color="green">194</FONT> final boolean forward = t &gt; t0;<a name="line.194"></a>
198 <FONT color="green">195</FONT> <a name="line.195"></a>
199 <FONT color="green">196</FONT> // initialize working arrays<a name="line.196"></a>
200 <FONT color="green">197</FONT> if (y != y0) {<a name="line.197"></a>
201 <FONT color="green">198</FONT> System.arraycopy(y0, 0, y, 0, n);<a name="line.198"></a>
202 <FONT color="green">199</FONT> }<a name="line.199"></a>
203 <FONT color="green">200</FONT> final double[] yDot = new double[n];<a name="line.200"></a>
204 <FONT color="green">201</FONT> final double[] yTmp = new double[y0.length];<a name="line.201"></a>
205 <FONT color="green">202</FONT> <a name="line.202"></a>
206 <FONT color="green">203</FONT> // set up an interpolator sharing the integrator arrays<a name="line.203"></a>
207 <FONT color="green">204</FONT> final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();<a name="line.204"></a>
208 <FONT color="green">205</FONT> interpolator.reinitialize(y, forward);<a name="line.205"></a>
209 <FONT color="green">206</FONT> final NordsieckStepInterpolator interpolatorTmp = new NordsieckStepInterpolator();<a name="line.206"></a>
210 <FONT color="green">207</FONT> interpolatorTmp.reinitialize(yTmp, forward);<a name="line.207"></a>
211 <FONT color="green">208</FONT> <a name="line.208"></a>
212 <FONT color="green">209</FONT> // set up integration control objects<a name="line.209"></a>
213 <FONT color="green">210</FONT> for (StepHandler handler : stepHandlers) {<a name="line.210"></a>
214 <FONT color="green">211</FONT> handler.reset();<a name="line.211"></a>
215 <FONT color="green">212</FONT> }<a name="line.212"></a>
216 <FONT color="green">213</FONT> CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);<a name="line.213"></a>
217 <FONT color="green">214</FONT> <a name="line.214"></a>
218 <FONT color="green">215</FONT> // compute the initial Nordsieck vector using the configured starter integrator<a name="line.215"></a>
219 <FONT color="green">216</FONT> start(t0, y, t);<a name="line.216"></a>
220 <FONT color="green">217</FONT> interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);<a name="line.217"></a>
221 <FONT color="green">218</FONT> interpolator.storeTime(stepStart);<a name="line.218"></a>
222 <FONT color="green">219</FONT> final int lastRow = nordsieck.getRowDimension() - 1;<a name="line.219"></a>
223 <FONT color="green">220</FONT> <a name="line.220"></a>
224 <FONT color="green">221</FONT> // reuse the step that was chosen by the starter integrator<a name="line.221"></a>
225 <FONT color="green">222</FONT> double hNew = stepSize;<a name="line.222"></a>
226 <FONT color="green">223</FONT> interpolator.rescale(hNew);<a name="line.223"></a>
227 <FONT color="green">224</FONT> <a name="line.224"></a>
228 <FONT color="green">225</FONT> boolean lastStep = false;<a name="line.225"></a>
229 <FONT color="green">226</FONT> while (!lastStep) {<a name="line.226"></a>
230 <FONT color="green">227</FONT> <a name="line.227"></a>
231 <FONT color="green">228</FONT> // shift all data<a name="line.228"></a>
232 <FONT color="green">229</FONT> interpolator.shift();<a name="line.229"></a>
233 <FONT color="green">230</FONT> <a name="line.230"></a>
234 <FONT color="green">231</FONT> double error = 0;<a name="line.231"></a>
235 <FONT color="green">232</FONT> for (boolean loop = true; loop;) {<a name="line.232"></a>
236 <FONT color="green">233</FONT> <a name="line.233"></a>
237 <FONT color="green">234</FONT> stepSize = hNew;<a name="line.234"></a>
238 <FONT color="green">235</FONT> <a name="line.235"></a>
239 <FONT color="green">236</FONT> // evaluate error using the last term of the Taylor expansion<a name="line.236"></a>
240 <FONT color="green">237</FONT> error = 0;<a name="line.237"></a>
241 <FONT color="green">238</FONT> for (int i = 0; i &lt; y0.length; ++i) {<a name="line.238"></a>
242 <FONT color="green">239</FONT> final double yScale = Math.abs(y[i]);<a name="line.239"></a>
243 <FONT color="green">240</FONT> final double tol = (vecAbsoluteTolerance == null) ?<a name="line.240"></a>
244 <FONT color="green">241</FONT> (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :<a name="line.241"></a>
245 <FONT color="green">242</FONT> (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);<a name="line.242"></a>
246 <FONT color="green">243</FONT> final double ratio = nordsieck.getEntry(lastRow, i) / tol;<a name="line.243"></a>
247 <FONT color="green">244</FONT> error += ratio * ratio;<a name="line.244"></a>
248 <FONT color="green">245</FONT> }<a name="line.245"></a>
249 <FONT color="green">246</FONT> error = Math.sqrt(error / y0.length);<a name="line.246"></a>
250 <FONT color="green">247</FONT> <a name="line.247"></a>
251 <FONT color="green">248</FONT> if (error &lt;= 1.0) {<a name="line.248"></a>
252 <FONT color="green">249</FONT> <a name="line.249"></a>
253 <FONT color="green">250</FONT> // predict a first estimate of the state at step end<a name="line.250"></a>
254 <FONT color="green">251</FONT> final double stepEnd = stepStart + stepSize;<a name="line.251"></a>
255 <FONT color="green">252</FONT> interpolator.setInterpolatedTime(stepEnd);<a name="line.252"></a>
256 <FONT color="green">253</FONT> System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length);<a name="line.253"></a>
257 <FONT color="green">254</FONT> <a name="line.254"></a>
258 <FONT color="green">255</FONT> // evaluate the derivative<a name="line.255"></a>
259 <FONT color="green">256</FONT> computeDerivatives(stepEnd, yTmp, yDot);<a name="line.256"></a>
260 <FONT color="green">257</FONT> <a name="line.257"></a>
261 <FONT color="green">258</FONT> // update Nordsieck vector<a name="line.258"></a>
262 <FONT color="green">259</FONT> final double[] predictedScaled = new double[y0.length];<a name="line.259"></a>
263 <FONT color="green">260</FONT> for (int j = 0; j &lt; y0.length; ++j) {<a name="line.260"></a>
264 <FONT color="green">261</FONT> predictedScaled[j] = stepSize * yDot[j];<a name="line.261"></a>
265 <FONT color="green">262</FONT> }<a name="line.262"></a>
266 <FONT color="green">263</FONT> final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);<a name="line.263"></a>
267 <FONT color="green">264</FONT> updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);<a name="line.264"></a>
268 <FONT color="green">265</FONT> <a name="line.265"></a>
269 <FONT color="green">266</FONT> // discrete events handling<a name="line.266"></a>
270 <FONT color="green">267</FONT> interpolatorTmp.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);<a name="line.267"></a>
271 <FONT color="green">268</FONT> interpolatorTmp.storeTime(stepStart);<a name="line.268"></a>
272 <FONT color="green">269</FONT> interpolatorTmp.shift();<a name="line.269"></a>
273 <FONT color="green">270</FONT> interpolatorTmp.storeTime(stepEnd);<a name="line.270"></a>
274 <FONT color="green">271</FONT> if (manager.evaluateStep(interpolatorTmp)) {<a name="line.271"></a>
275 <FONT color="green">272</FONT> final double dt = manager.getEventTime() - stepStart;<a name="line.272"></a>
276 <FONT color="green">273</FONT> if (Math.abs(dt) &lt;= Math.ulp(stepStart)) {<a name="line.273"></a>
277 <FONT color="green">274</FONT> // we cannot simply truncate the step, reject the current computation<a name="line.274"></a>
278 <FONT color="green">275</FONT> // and let the loop compute another state with the truncated step.<a name="line.275"></a>
279 <FONT color="green">276</FONT> // it is so small (much probably exactly 0 due to limited accuracy)<a name="line.276"></a>
280 <FONT color="green">277</FONT> // that the code above would fail handling it.<a name="line.277"></a>
281 <FONT color="green">278</FONT> // So we set up an artificial 0 size step by copying states<a name="line.278"></a>
282 <FONT color="green">279</FONT> interpolator.storeTime(stepStart);<a name="line.279"></a>
283 <FONT color="green">280</FONT> System.arraycopy(y, 0, yTmp, 0, y0.length);<a name="line.280"></a>
284 <FONT color="green">281</FONT> hNew = 0;<a name="line.281"></a>
285 <FONT color="green">282</FONT> stepSize = 0;<a name="line.282"></a>
286 <FONT color="green">283</FONT> loop = false;<a name="line.283"></a>
287 <FONT color="green">284</FONT> } else {<a name="line.284"></a>
288 <FONT color="green">285</FONT> // reject the step to match exactly the next switch time<a name="line.285"></a>
289 <FONT color="green">286</FONT> hNew = dt;<a name="line.286"></a>
290 <FONT color="green">287</FONT> interpolator.rescale(hNew);<a name="line.287"></a>
291 <FONT color="green">288</FONT> }<a name="line.288"></a>
292 <FONT color="green">289</FONT> } else {<a name="line.289"></a>
293 <FONT color="green">290</FONT> // accept the step<a name="line.290"></a>
294 <FONT color="green">291</FONT> scaled = predictedScaled;<a name="line.291"></a>
295 <FONT color="green">292</FONT> nordsieck = nordsieckTmp;<a name="line.292"></a>
296 <FONT color="green">293</FONT> interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);<a name="line.293"></a>
297 <FONT color="green">294</FONT> loop = false;<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> } else {<a name="line.297"></a>
301 <FONT color="green">298</FONT> // reject the step and attempt to reduce error by stepsize control<a name="line.298"></a>
302 <FONT color="green">299</FONT> final double factor = computeStepGrowShrinkFactor(error);<a name="line.299"></a>
303 <FONT color="green">300</FONT> hNew = filterStep(stepSize * factor, forward, false);<a name="line.300"></a>
304 <FONT color="green">301</FONT> interpolator.rescale(hNew);<a name="line.301"></a>
305 <FONT color="green">302</FONT> }<a name="line.302"></a>
306 <FONT color="green">303</FONT> <a name="line.303"></a>
307 <FONT color="green">304</FONT> }<a name="line.304"></a>
308 <FONT color="green">305</FONT> <a name="line.305"></a>
309 <FONT color="green">306</FONT> // the step has been accepted (may have been truncated)<a name="line.306"></a>
310 <FONT color="green">307</FONT> final double nextStep = stepStart + stepSize;<a name="line.307"></a>
311 <FONT color="green">308</FONT> System.arraycopy(yTmp, 0, y, 0, n);<a name="line.308"></a>
312 <FONT color="green">309</FONT> interpolator.storeTime(nextStep);<a name="line.309"></a>
313 <FONT color="green">310</FONT> manager.stepAccepted(nextStep, y);<a name="line.310"></a>
314 <FONT color="green">311</FONT> lastStep = manager.stop();<a name="line.311"></a>
315 <FONT color="green">312</FONT> <a name="line.312"></a>
316 <FONT color="green">313</FONT> // provide the step data to the step handler<a name="line.313"></a>
317 <FONT color="green">314</FONT> for (StepHandler handler : stepHandlers) {<a name="line.314"></a>
318 <FONT color="green">315</FONT> interpolator.setInterpolatedTime(nextStep);<a name="line.315"></a>
319 <FONT color="green">316</FONT> handler.handleStep(interpolator, lastStep);<a name="line.316"></a>
320 <FONT color="green">317</FONT> }<a name="line.317"></a>
321 <FONT color="green">318</FONT> stepStart = nextStep;<a name="line.318"></a>
322 <FONT color="green">319</FONT> <a name="line.319"></a>
323 <FONT color="green">320</FONT> if (!lastStep &amp;&amp; manager.reset(stepStart, y)) {<a name="line.320"></a>
324 <FONT color="green">321</FONT> <a name="line.321"></a>
325 <FONT color="green">322</FONT> // some events handler has triggered changes that<a name="line.322"></a>
326 <FONT color="green">323</FONT> // invalidate the derivatives, we need to restart from scratch<a name="line.323"></a>
327 <FONT color="green">324</FONT> start(stepStart, y, t);<a name="line.324"></a>
328 <FONT color="green">325</FONT> interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);<a name="line.325"></a>
329 <FONT color="green">326</FONT> <a name="line.326"></a>
330 <FONT color="green">327</FONT> }<a name="line.327"></a>
331 <FONT color="green">328</FONT> <a name="line.328"></a>
332 <FONT color="green">329</FONT> if (! lastStep) {<a name="line.329"></a>
333 <FONT color="green">330</FONT> // in some rare cases we may get here with stepSize = 0, for example<a name="line.330"></a>
334 <FONT color="green">331</FONT> // when an event occurs at integration start, reducing the first step<a name="line.331"></a>
335 <FONT color="green">332</FONT> // to zero; we have to reset the step to some safe non zero value<a name="line.332"></a>
336 <FONT color="green">333</FONT> stepSize = filterStep(stepSize, forward, true);<a name="line.333"></a>
337 <FONT color="green">334</FONT> <a name="line.334"></a>
338 <FONT color="green">335</FONT> // stepsize control for next step<a name="line.335"></a>
339 <FONT color="green">336</FONT> final double factor = computeStepGrowShrinkFactor(error);<a name="line.336"></a>
340 <FONT color="green">337</FONT> final double scaledH = stepSize * factor;<a name="line.337"></a>
341 <FONT color="green">338</FONT> final double nextT = stepStart + scaledH;<a name="line.338"></a>
342 <FONT color="green">339</FONT> final boolean nextIsLast = forward ? (nextT &gt;= t) : (nextT &lt;= t);<a name="line.339"></a>
343 <FONT color="green">340</FONT> hNew = filterStep(scaledH, forward, nextIsLast);<a name="line.340"></a>
344 <FONT color="green">341</FONT> interpolator.rescale(hNew);<a name="line.341"></a>
345 <FONT color="green">342</FONT> }<a name="line.342"></a>
346 <FONT color="green">343</FONT> <a name="line.343"></a>
347 <FONT color="green">344</FONT> }<a name="line.344"></a>
348 <FONT color="green">345</FONT> <a name="line.345"></a>
349 <FONT color="green">346</FONT> final double stopTime = stepStart;<a name="line.346"></a>
350 <FONT color="green">347</FONT> stepStart = Double.NaN;<a name="line.347"></a>
351 <FONT color="green">348</FONT> stepSize = Double.NaN;<a name="line.348"></a>
352 <FONT color="green">349</FONT> return stopTime;<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> }<a name="line.353"></a>
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
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417 </PRE>
418 </BODY>
419 </HTML>