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

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.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 java.util.Arrays;<a name="line.20"></a>
24 <FONT color="green">021</FONT> import java.util.HashMap;<a name="line.21"></a>
25 <FONT color="green">022</FONT> import java.util.Map;<a name="line.22"></a>
26 <FONT color="green">023</FONT> <a name="line.23"></a>
27 <FONT color="green">024</FONT> import org.apache.commons.math.fraction.BigFraction;<a name="line.24"></a>
28 <FONT color="green">025</FONT> import org.apache.commons.math.linear.Array2DRowFieldMatrix;<a name="line.25"></a>
29 <FONT color="green">026</FONT> import org.apache.commons.math.linear.Array2DRowRealMatrix;<a name="line.26"></a>
30 <FONT color="green">027</FONT> import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor;<a name="line.27"></a>
31 <FONT color="green">028</FONT> import org.apache.commons.math.linear.FieldDecompositionSolver;<a name="line.28"></a>
32 <FONT color="green">029</FONT> import org.apache.commons.math.linear.FieldLUDecompositionImpl;<a name="line.29"></a>
33 <FONT color="green">030</FONT> import org.apache.commons.math.linear.FieldMatrix;<a name="line.30"></a>
34 <FONT color="green">031</FONT> import org.apache.commons.math.linear.MatrixUtils;<a name="line.31"></a>
35 <FONT color="green">032</FONT> <a name="line.32"></a>
36 <FONT color="green">033</FONT> /** Transformer to Nordsieck vectors for Adams integrators.<a name="line.33"></a>
37 <FONT color="green">034</FONT> * &lt;p&gt;This class i used by {@link AdamsBashforthIntegrator Adams-Bashforth} and<a name="line.34"></a>
38 <FONT color="green">035</FONT> * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between<a name="line.35"></a>
39 <FONT color="green">036</FONT> * classical representation with several previous first derivatives and Nordsieck<a name="line.36"></a>
40 <FONT color="green">037</FONT> * representation with higher order scaled derivatives.&lt;/p&gt;<a name="line.37"></a>
41 <FONT color="green">038</FONT> *<a name="line.38"></a>
42 <FONT color="green">039</FONT> * &lt;p&gt;We define scaled derivatives s&lt;sub&gt;i&lt;/sub&gt;(n) at step n as:<a name="line.39"></a>
43 <FONT color="green">040</FONT> * &lt;pre&gt;<a name="line.40"></a>
44 <FONT color="green">041</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.41"></a>
45 <FONT color="green">042</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.42"></a>
46 <FONT color="green">043</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.43"></a>
47 <FONT color="green">044</FONT> * ...<a name="line.44"></a>
48 <FONT color="green">045</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.45"></a>
49 <FONT color="green">046</FONT> * &lt;/pre&gt;&lt;/p&gt;<a name="line.46"></a>
50 <FONT color="green">047</FONT> *<a name="line.47"></a>
51 <FONT color="green">048</FONT> * &lt;p&gt;With the previous definition, the classical representation of multistep methods<a name="line.48"></a>
52 <FONT color="green">049</FONT> * uses first derivatives only, i.e. it handles y&lt;sub&gt;n&lt;/sub&gt;, s&lt;sub&gt;1&lt;/sub&gt;(n) and<a name="line.49"></a>
53 <FONT color="green">050</FONT> * q&lt;sub&gt;n&lt;/sub&gt; where q&lt;sub&gt;n&lt;/sub&gt; is defined as:<a name="line.50"></a>
54 <FONT color="green">051</FONT> * &lt;pre&gt;<a name="line.51"></a>
55 <FONT color="green">052</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.52"></a>
56 <FONT color="green">053</FONT> * &lt;/pre&gt;<a name="line.53"></a>
57 <FONT color="green">054</FONT> * (we omit the k index in the notation for clarity).&lt;/p&gt;<a name="line.54"></a>
58 <FONT color="green">055</FONT> *<a name="line.55"></a>
59 <FONT color="green">056</FONT> * &lt;p&gt;Another possible representation uses the Nordsieck vector with<a name="line.56"></a>
60 <FONT color="green">057</FONT> * higher degrees scaled derivatives all taken at the same step, i.e it handles y&lt;sub&gt;n&lt;/sub&gt;,<a name="line.57"></a>
61 <FONT color="green">058</FONT> * 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.58"></a>
62 <FONT color="green">059</FONT> * &lt;pre&gt;<a name="line.59"></a>
63 <FONT color="green">060</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.60"></a>
64 <FONT color="green">061</FONT> * &lt;/pre&gt;<a name="line.61"></a>
65 <FONT color="green">062</FONT> * (here again we omit the k index in the notation for clarity)<a name="line.62"></a>
66 <FONT color="green">063</FONT> * &lt;/p&gt;<a name="line.63"></a>
67 <FONT color="green">064</FONT> *<a name="line.64"></a>
68 <FONT color="green">065</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.65"></a>
69 <FONT color="green">066</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.66"></a>
70 <FONT color="green">067</FONT> * for degree k polynomials.<a name="line.67"></a>
71 <FONT color="green">068</FONT> * &lt;pre&gt;<a name="line.68"></a>
72 <FONT color="green">069</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.69"></a>
73 <FONT color="green">070</FONT> * &lt;/pre&gt;<a name="line.70"></a>
74 <FONT color="green">071</FONT> * The previous formula can be used with several values for i to compute the transform between<a name="line.71"></a>
75 <FONT color="green">072</FONT> * classical representation and Nordsieck vector at step end. The transform between r&lt;sub&gt;n&lt;/sub&gt;<a name="line.72"></a>
76 <FONT color="green">073</FONT> * and q&lt;sub&gt;n&lt;/sub&gt; resulting from the Taylor series formulas above is:<a name="line.73"></a>
77 <FONT color="green">074</FONT> * &lt;pre&gt;<a name="line.74"></a>
78 <FONT color="green">075</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.75"></a>
79 <FONT color="green">076</FONT> * &lt;/pre&gt;<a name="line.76"></a>
80 <FONT color="green">077</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.77"></a>
81 <FONT color="green">078</FONT> * with the j (-i)&lt;sup&gt;j-1&lt;/sup&gt; terms:<a name="line.78"></a>
82 <FONT color="green">079</FONT> * &lt;pre&gt;<a name="line.79"></a>
83 <FONT color="green">080</FONT> * [ -2 3 -4 5 ... ]<a name="line.80"></a>
84 <FONT color="green">081</FONT> * [ -4 12 -32 80 ... ]<a name="line.81"></a>
85 <FONT color="green">082</FONT> * P = [ -6 27 -108 405 ... ]<a name="line.82"></a>
86 <FONT color="green">083</FONT> * [ -8 48 -256 1280 ... ]<a name="line.83"></a>
87 <FONT color="green">084</FONT> * [ ... ]<a name="line.84"></a>
88 <FONT color="green">085</FONT> * &lt;/pre&gt;&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;Changing -i into +i in the formula above can be used to compute a similar transform between<a name="line.87"></a>
91 <FONT color="green">088</FONT> * classical representation and Nordsieck vector at step start. The resulting matrix is simply<a name="line.88"></a>
92 <FONT color="green">089</FONT> * the absolute value of matrix P.&lt;/p&gt;<a name="line.89"></a>
93 <FONT color="green">090</FONT> *<a name="line.90"></a>
94 <FONT color="green">091</FONT> * &lt;p&gt;For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector<a name="line.91"></a>
95 <FONT color="green">092</FONT> * at step n+1 is computed from the Nordsieck vector at step n as follows:<a name="line.92"></a>
96 <FONT color="green">093</FONT> * &lt;ul&gt;<a name="line.93"></a>
97 <FONT color="green">094</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.94"></a>
98 <FONT color="green">095</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.95"></a>
99 <FONT color="green">096</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.96"></a>
100 <FONT color="green">097</FONT> * &lt;/ul&gt;<a name="line.97"></a>
101 <FONT color="green">098</FONT> * where A is a rows shifting matrix (the lower left part is an identity matrix):<a name="line.98"></a>
102 <FONT color="green">099</FONT> * &lt;pre&gt;<a name="line.99"></a>
103 <FONT color="green">100</FONT> * [ 0 0 ... 0 0 | 0 ]<a name="line.100"></a>
104 <FONT color="green">101</FONT> * [ ---------------+---]<a name="line.101"></a>
105 <FONT color="green">102</FONT> * [ 1 0 ... 0 0 | 0 ]<a name="line.102"></a>
106 <FONT color="green">103</FONT> * A = [ 0 1 ... 0 0 | 0 ]<a name="line.103"></a>
107 <FONT color="green">104</FONT> * [ ... | 0 ]<a name="line.104"></a>
108 <FONT color="green">105</FONT> * [ 0 0 ... 1 0 | 0 ]<a name="line.105"></a>
109 <FONT color="green">106</FONT> * [ 0 0 ... 0 1 | 0 ]<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;For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector<a name="line.109"></a>
113 <FONT color="green">110</FONT> * at step n+1 is computed from the Nordsieck vector at step n as follows:<a name="line.110"></a>
114 <FONT color="green">111</FONT> * &lt;ul&gt;<a name="line.111"></a>
115 <FONT color="green">112</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.112"></a>
116 <FONT color="green">113</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.113"></a>
117 <FONT color="green">114</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.114"></a>
118 <FONT color="green">115</FONT> * &lt;/ul&gt;<a name="line.115"></a>
119 <FONT color="green">116</FONT> * From this predicted vector, the corrected vector is computed as follows:<a name="line.116"></a>
120 <FONT color="green">117</FONT> * &lt;ul&gt;<a name="line.117"></a>
121 <FONT color="green">118</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+1) + [ -1 +1 -1 +1 ... &amp;plusmn;1 ] r&lt;sub&gt;n+1&lt;/sub&gt;&lt;/li&gt;<a name="line.118"></a>
122 <FONT color="green">119</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.119"></a>
123 <FONT color="green">120</FONT> * &lt;li&gt;r&lt;sub&gt;n+1&lt;/sub&gt; = R&lt;sub&gt;n+1&lt;/sub&gt; + (s&lt;sub&gt;1&lt;/sub&gt;(n+1) - S&lt;sub&gt;1&lt;/sub&gt;(n+1)) P&lt;sup&gt;-1&lt;/sup&gt; u&lt;/li&gt;<a name="line.120"></a>
124 <FONT color="green">121</FONT> * &lt;/ul&gt;<a name="line.121"></a>
125 <FONT color="green">122</FONT> * where the upper case Y&lt;sub&gt;n+1&lt;/sub&gt;, S&lt;sub&gt;1&lt;/sub&gt;(n+1) and R&lt;sub&gt;n+1&lt;/sub&gt; represent the<a name="line.122"></a>
126 <FONT color="green">123</FONT> * predicted states whereas the lower case y&lt;sub&gt;n+1&lt;/sub&gt;, s&lt;sub&gt;n+1&lt;/sub&gt; and r&lt;sub&gt;n+1&lt;/sub&gt;<a name="line.123"></a>
127 <FONT color="green">124</FONT> * represent the corrected states.&lt;/p&gt;<a name="line.124"></a>
128 <FONT color="green">125</FONT> *<a name="line.125"></a>
129 <FONT color="green">126</FONT> * &lt;p&gt;We observe that both methods use similar update formulas. In both cases a P&lt;sup&gt;-1&lt;/sup&gt;u<a name="line.126"></a>
130 <FONT color="green">127</FONT> * vector and a P&lt;sup&gt;-1&lt;/sup&gt; A P matrix are used that do not depend on the state,<a name="line.127"></a>
131 <FONT color="green">128</FONT> * they only depend on k. This class handles these transformations.&lt;/p&gt;<a name="line.128"></a>
132 <FONT color="green">129</FONT> *<a name="line.129"></a>
133 <FONT color="green">130</FONT> * @version $Revision: 810196 $ $Date: 2009-09-01 15:47:46 -0400 (Tue, 01 Sep 2009) $<a name="line.130"></a>
134 <FONT color="green">131</FONT> * @since 2.0<a name="line.131"></a>
135 <FONT color="green">132</FONT> */<a name="line.132"></a>
136 <FONT color="green">133</FONT> public class AdamsNordsieckTransformer {<a name="line.133"></a>
137 <FONT color="green">134</FONT> <a name="line.134"></a>
138 <FONT color="green">135</FONT> /** Cache for already computed coefficients. */<a name="line.135"></a>
139 <FONT color="green">136</FONT> private static final Map&lt;Integer, AdamsNordsieckTransformer&gt; CACHE =<a name="line.136"></a>
140 <FONT color="green">137</FONT> new HashMap&lt;Integer, AdamsNordsieckTransformer&gt;();<a name="line.137"></a>
141 <FONT color="green">138</FONT> <a name="line.138"></a>
142 <FONT color="green">139</FONT> /** Initialization matrix for the higher order derivatives wrt y'', y''' ... */<a name="line.139"></a>
143 <FONT color="green">140</FONT> private final Array2DRowRealMatrix initialization;<a name="line.140"></a>
144 <FONT color="green">141</FONT> <a name="line.141"></a>
145 <FONT color="green">142</FONT> /** Update matrix for the higher order derivatives h&lt;sup&gt;2&lt;/sup&gt;/2y'', h&lt;sup&gt;3&lt;/sup&gt;/6 y''' ... */<a name="line.142"></a>
146 <FONT color="green">143</FONT> private final Array2DRowRealMatrix update;<a name="line.143"></a>
147 <FONT color="green">144</FONT> <a name="line.144"></a>
148 <FONT color="green">145</FONT> /** Update coefficients of the higher order derivatives wrt y'. */<a name="line.145"></a>
149 <FONT color="green">146</FONT> private final double[] c1;<a name="line.146"></a>
150 <FONT color="green">147</FONT> <a name="line.147"></a>
151 <FONT color="green">148</FONT> /** Simple constructor.<a name="line.148"></a>
152 <FONT color="green">149</FONT> * @param nSteps number of steps of the multistep method<a name="line.149"></a>
153 <FONT color="green">150</FONT> * (excluding the one being computed)<a name="line.150"></a>
154 <FONT color="green">151</FONT> */<a name="line.151"></a>
155 <FONT color="green">152</FONT> private AdamsNordsieckTransformer(final int nSteps) {<a name="line.152"></a>
156 <FONT color="green">153</FONT> <a name="line.153"></a>
157 <FONT color="green">154</FONT> // compute exact coefficients<a name="line.154"></a>
158 <FONT color="green">155</FONT> FieldMatrix&lt;BigFraction&gt; bigP = buildP(nSteps);<a name="line.155"></a>
159 <FONT color="green">156</FONT> FieldDecompositionSolver&lt;BigFraction&gt; pSolver =<a name="line.156"></a>
160 <FONT color="green">157</FONT> new FieldLUDecompositionImpl&lt;BigFraction&gt;(bigP).getSolver();<a name="line.157"></a>
161 <FONT color="green">158</FONT> <a name="line.158"></a>
162 <FONT color="green">159</FONT> BigFraction[] u = new BigFraction[nSteps];<a name="line.159"></a>
163 <FONT color="green">160</FONT> Arrays.fill(u, BigFraction.ONE);<a name="line.160"></a>
164 <FONT color="green">161</FONT> BigFraction[] bigC1 = pSolver.solve(u);<a name="line.161"></a>
165 <FONT color="green">162</FONT> <a name="line.162"></a>
166 <FONT color="green">163</FONT> // update coefficients are computed by combining transform from<a name="line.163"></a>
167 <FONT color="green">164</FONT> // Nordsieck to multistep, then shifting rows to represent step advance<a name="line.164"></a>
168 <FONT color="green">165</FONT> // then applying inverse transform<a name="line.165"></a>
169 <FONT color="green">166</FONT> BigFraction[][] shiftedP = bigP.getData();<a name="line.166"></a>
170 <FONT color="green">167</FONT> for (int i = shiftedP.length - 1; i &gt; 0; --i) {<a name="line.167"></a>
171 <FONT color="green">168</FONT> // shift rows<a name="line.168"></a>
172 <FONT color="green">169</FONT> shiftedP[i] = shiftedP[i - 1];<a name="line.169"></a>
173 <FONT color="green">170</FONT> }<a name="line.170"></a>
174 <FONT color="green">171</FONT> shiftedP[0] = new BigFraction[nSteps];<a name="line.171"></a>
175 <FONT color="green">172</FONT> Arrays.fill(shiftedP[0], BigFraction.ZERO);<a name="line.172"></a>
176 <FONT color="green">173</FONT> FieldMatrix&lt;BigFraction&gt; bigMSupdate =<a name="line.173"></a>
177 <FONT color="green">174</FONT> pSolver.solve(new Array2DRowFieldMatrix&lt;BigFraction&gt;(shiftedP, false));<a name="line.174"></a>
178 <FONT color="green">175</FONT> <a name="line.175"></a>
179 <FONT color="green">176</FONT> // initialization coefficients, computed from a R matrix = abs(P)<a name="line.176"></a>
180 <FONT color="green">177</FONT> bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor&lt;BigFraction&gt;(BigFraction.ZERO) {<a name="line.177"></a>
181 <FONT color="green">178</FONT> /** {@inheritDoc} */<a name="line.178"></a>
182 <FONT color="green">179</FONT> @Override<a name="line.179"></a>
183 <FONT color="green">180</FONT> public BigFraction visit(int row, int column, BigFraction value) {<a name="line.180"></a>
184 <FONT color="green">181</FONT> return ((column &amp; 0x1) == 0x1) ? value : value.negate();<a name="line.181"></a>
185 <FONT color="green">182</FONT> }<a name="line.182"></a>
186 <FONT color="green">183</FONT> });<a name="line.183"></a>
187 <FONT color="green">184</FONT> FieldMatrix&lt;BigFraction&gt; bigRInverse =<a name="line.184"></a>
188 <FONT color="green">185</FONT> new FieldLUDecompositionImpl&lt;BigFraction&gt;(bigP).getSolver().getInverse();<a name="line.185"></a>
189 <FONT color="green">186</FONT> <a name="line.186"></a>
190 <FONT color="green">187</FONT> // convert coefficients to double<a name="line.187"></a>
191 <FONT color="green">188</FONT> initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigRInverse);<a name="line.188"></a>
192 <FONT color="green">189</FONT> update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);<a name="line.189"></a>
193 <FONT color="green">190</FONT> c1 = new double[nSteps];<a name="line.190"></a>
194 <FONT color="green">191</FONT> for (int i = 0; i &lt; nSteps; ++i) {<a name="line.191"></a>
195 <FONT color="green">192</FONT> c1[i] = bigC1[i].doubleValue();<a name="line.192"></a>
196 <FONT color="green">193</FONT> }<a name="line.193"></a>
197 <FONT color="green">194</FONT> <a name="line.194"></a>
198 <FONT color="green">195</FONT> }<a name="line.195"></a>
199 <FONT color="green">196</FONT> <a name="line.196"></a>
200 <FONT color="green">197</FONT> /** Get the Nordsieck transformer for a given number of steps.<a name="line.197"></a>
201 <FONT color="green">198</FONT> * @param nSteps number of steps of the multistep method<a name="line.198"></a>
202 <FONT color="green">199</FONT> * (excluding the one being computed)<a name="line.199"></a>
203 <FONT color="green">200</FONT> * @return Nordsieck transformer for the specified number of steps<a name="line.200"></a>
204 <FONT color="green">201</FONT> */<a name="line.201"></a>
205 <FONT color="green">202</FONT> public static AdamsNordsieckTransformer getInstance(final int nSteps) {<a name="line.202"></a>
206 <FONT color="green">203</FONT> synchronized(CACHE) {<a name="line.203"></a>
207 <FONT color="green">204</FONT> AdamsNordsieckTransformer t = CACHE.get(nSteps);<a name="line.204"></a>
208 <FONT color="green">205</FONT> if (t == null) {<a name="line.205"></a>
209 <FONT color="green">206</FONT> t = new AdamsNordsieckTransformer(nSteps);<a name="line.206"></a>
210 <FONT color="green">207</FONT> CACHE.put(nSteps, t);<a name="line.207"></a>
211 <FONT color="green">208</FONT> }<a name="line.208"></a>
212 <FONT color="green">209</FONT> return t;<a name="line.209"></a>
213 <FONT color="green">210</FONT> }<a name="line.210"></a>
214 <FONT color="green">211</FONT> }<a name="line.211"></a>
215 <FONT color="green">212</FONT> <a name="line.212"></a>
216 <FONT color="green">213</FONT> /** Get the number of steps of the method<a name="line.213"></a>
217 <FONT color="green">214</FONT> * (excluding the one being computed).<a name="line.214"></a>
218 <FONT color="green">215</FONT> * @return number of steps of the method<a name="line.215"></a>
219 <FONT color="green">216</FONT> * (excluding the one being computed)<a name="line.216"></a>
220 <FONT color="green">217</FONT> */<a name="line.217"></a>
221 <FONT color="green">218</FONT> public int getNSteps() {<a name="line.218"></a>
222 <FONT color="green">219</FONT> return c1.length;<a name="line.219"></a>
223 <FONT color="green">220</FONT> }<a name="line.220"></a>
224 <FONT color="green">221</FONT> <a name="line.221"></a>
225 <FONT color="green">222</FONT> /** Build the P matrix.<a name="line.222"></a>
226 <FONT color="green">223</FONT> * &lt;p&gt;The P matrix general terms are shifted j (-i)&lt;sup&gt;j-1&lt;/sup&gt; terms:<a name="line.223"></a>
227 <FONT color="green">224</FONT> * &lt;pre&gt;<a name="line.224"></a>
228 <FONT color="green">225</FONT> * [ -2 3 -4 5 ... ]<a name="line.225"></a>
229 <FONT color="green">226</FONT> * [ -4 12 -32 80 ... ]<a name="line.226"></a>
230 <FONT color="green">227</FONT> * P = [ -6 27 -108 405 ... ]<a name="line.227"></a>
231 <FONT color="green">228</FONT> * [ -8 48 -256 1280 ... ]<a name="line.228"></a>
232 <FONT color="green">229</FONT> * [ ... ]<a name="line.229"></a>
233 <FONT color="green">230</FONT> * &lt;/pre&gt;&lt;/p&gt;<a name="line.230"></a>
234 <FONT color="green">231</FONT> * @param nSteps number of steps of the multistep method<a name="line.231"></a>
235 <FONT color="green">232</FONT> * (excluding the one being computed)<a name="line.232"></a>
236 <FONT color="green">233</FONT> * @return P matrix<a name="line.233"></a>
237 <FONT color="green">234</FONT> */<a name="line.234"></a>
238 <FONT color="green">235</FONT> private FieldMatrix&lt;BigFraction&gt; buildP(final int nSteps) {<a name="line.235"></a>
239 <FONT color="green">236</FONT> <a name="line.236"></a>
240 <FONT color="green">237</FONT> final BigFraction[][] pData = new BigFraction[nSteps][nSteps];<a name="line.237"></a>
241 <FONT color="green">238</FONT> <a name="line.238"></a>
242 <FONT color="green">239</FONT> for (int i = 0; i &lt; pData.length; ++i) {<a name="line.239"></a>
243 <FONT color="green">240</FONT> // build the P matrix elements from Taylor series formulas<a name="line.240"></a>
244 <FONT color="green">241</FONT> final BigFraction[] pI = pData[i];<a name="line.241"></a>
245 <FONT color="green">242</FONT> final int factor = -(i + 1);<a name="line.242"></a>
246 <FONT color="green">243</FONT> int aj = factor;<a name="line.243"></a>
247 <FONT color="green">244</FONT> for (int j = 0; j &lt; pI.length; ++j) {<a name="line.244"></a>
248 <FONT color="green">245</FONT> pI[j] = new BigFraction(aj * (j + 2));<a name="line.245"></a>
249 <FONT color="green">246</FONT> aj *= factor;<a name="line.246"></a>
250 <FONT color="green">247</FONT> }<a name="line.247"></a>
251 <FONT color="green">248</FONT> }<a name="line.248"></a>
252 <FONT color="green">249</FONT> <a name="line.249"></a>
253 <FONT color="green">250</FONT> return new Array2DRowFieldMatrix&lt;BigFraction&gt;(pData, false);<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> <a name="line.253"></a>
257 <FONT color="green">254</FONT> /** Initialize the high order scaled derivatives at step start.<a name="line.254"></a>
258 <FONT color="green">255</FONT> * @param first first scaled derivative at step start<a name="line.255"></a>
259 <FONT color="green">256</FONT> * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)<a name="line.256"></a>
260 <FONT color="green">257</FONT> * will be modified<a name="line.257"></a>
261 <FONT color="green">258</FONT> * @return high order derivatives at step start<a name="line.258"></a>
262 <FONT color="green">259</FONT> */<a name="line.259"></a>
263 <FONT color="green">260</FONT> public Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,<a name="line.260"></a>
264 <FONT color="green">261</FONT> final double[][] multistep) {<a name="line.261"></a>
265 <FONT color="green">262</FONT> for (int i = 0; i &lt; multistep.length; ++i) {<a name="line.262"></a>
266 <FONT color="green">263</FONT> final double[] msI = multistep[i];<a name="line.263"></a>
267 <FONT color="green">264</FONT> for (int j = 0; j &lt; first.length; ++j) {<a name="line.264"></a>
268 <FONT color="green">265</FONT> msI[j] -= first[j];<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> return initialization.multiply(new Array2DRowRealMatrix(multistep, false));<a name="line.268"></a>
272 <FONT color="green">269</FONT> }<a name="line.269"></a>
273 <FONT color="green">270</FONT> <a name="line.270"></a>
274 <FONT color="green">271</FONT> /** Update the high order scaled derivatives for Adams integrators (phase 1).<a name="line.271"></a>
275 <FONT color="green">272</FONT> * &lt;p&gt;The complete update of high order derivatives has a form similar to:<a name="line.272"></a>
276 <FONT color="green">273</FONT> * &lt;pre&gt;<a name="line.273"></a>
277 <FONT color="green">274</FONT> * 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;<a name="line.274"></a>
278 <FONT color="green">275</FONT> * &lt;/pre&gt;<a name="line.275"></a>
279 <FONT color="green">276</FONT> * this method computes the P&lt;sup&gt;-1&lt;/sup&gt; A P r&lt;sub&gt;n&lt;/sub&gt; part.&lt;/p&gt;<a name="line.276"></a>
280 <FONT color="green">277</FONT> * @param highOrder high order scaled derivatives<a name="line.277"></a>
281 <FONT color="green">278</FONT> * (h&lt;sup&gt;2&lt;/sup&gt;/2 y'', ... h&lt;sup&gt;k&lt;/sup&gt;/k! y(k))<a name="line.278"></a>
282 <FONT color="green">279</FONT> * @return updated high order derivatives<a name="line.279"></a>
283 <FONT color="green">280</FONT> * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)<a name="line.280"></a>
284 <FONT color="green">281</FONT> */<a name="line.281"></a>
285 <FONT color="green">282</FONT> public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {<a name="line.282"></a>
286 <FONT color="green">283</FONT> return update.multiply(highOrder);<a name="line.283"></a>
287 <FONT color="green">284</FONT> }<a name="line.284"></a>
288 <FONT color="green">285</FONT> <a name="line.285"></a>
289 <FONT color="green">286</FONT> /** Update the high order scaled derivatives Adams integrators (phase 2).<a name="line.286"></a>
290 <FONT color="green">287</FONT> * &lt;p&gt;The complete update of high order derivatives has a form similar to:<a name="line.287"></a>
291 <FONT color="green">288</FONT> * &lt;pre&gt;<a name="line.288"></a>
292 <FONT color="green">289</FONT> * 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;<a name="line.289"></a>
293 <FONT color="green">290</FONT> * &lt;/pre&gt;<a name="line.290"></a>
294 <FONT color="green">291</FONT> * this method computes the (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 part.&lt;/p&gt;<a name="line.291"></a>
295 <FONT color="green">292</FONT> * &lt;p&gt;Phase 1 of the update must already have been performed.&lt;/p&gt;<a name="line.292"></a>
296 <FONT color="green">293</FONT> * @param start first order scaled derivatives at step start<a name="line.293"></a>
297 <FONT color="green">294</FONT> * @param end first order scaled derivatives at step end<a name="line.294"></a>
298 <FONT color="green">295</FONT> * @param highOrder high order scaled derivatives, will be modified<a name="line.295"></a>
299 <FONT color="green">296</FONT> * (h&lt;sup&gt;2&lt;/sup&gt;/2 y'', ... h&lt;sup&gt;k&lt;/sup&gt;/k! y(k))<a name="line.296"></a>
300 <FONT color="green">297</FONT> * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)<a name="line.297"></a>
301 <FONT color="green">298</FONT> */<a name="line.298"></a>
302 <FONT color="green">299</FONT> public void updateHighOrderDerivativesPhase2(final double[] start,<a name="line.299"></a>
303 <FONT color="green">300</FONT> final double[] end,<a name="line.300"></a>
304 <FONT color="green">301</FONT> final Array2DRowRealMatrix highOrder) {<a name="line.301"></a>
305 <FONT color="green">302</FONT> final double[][] data = highOrder.getDataRef();<a name="line.302"></a>
306 <FONT color="green">303</FONT> for (int i = 0; i &lt; data.length; ++i) {<a name="line.303"></a>
307 <FONT color="green">304</FONT> final double[] dataI = data[i];<a name="line.304"></a>
308 <FONT color="green">305</FONT> final double c1I = c1[i];<a name="line.305"></a>
309 <FONT color="green">306</FONT> for (int j = 0; j &lt; dataI.length; ++j) {<a name="line.306"></a>
310 <FONT color="green">307</FONT> dataI[j] += c1I * (start[j] - end[j]);<a name="line.307"></a>
311 <FONT color="green">308</FONT> }<a name="line.308"></a>
312 <FONT color="green">309</FONT> }<a name="line.309"></a>
313 <FONT color="green">310</FONT> }<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
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
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 </PRE>
377 </BODY>
378 </HTML>