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

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.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> <a name="line.20"></a>
24 <FONT color="green">021</FONT> /**<a name="line.21"></a>
25 <FONT color="green">022</FONT> * This class implements the 8(5,3) Dormand-Prince integrator for Ordinary<a name="line.22"></a>
26 <FONT color="green">023</FONT> * Differential Equations.<a name="line.23"></a>
27 <FONT color="green">024</FONT> *<a name="line.24"></a>
28 <FONT color="green">025</FONT> * &lt;p&gt;This integrator is an embedded Runge-Kutta integrator<a name="line.25"></a>
29 <FONT color="green">026</FONT> * of order 8(5,3) used in local extrapolation mode (i.e. the solution<a name="line.26"></a>
30 <FONT color="green">027</FONT> * is computed using the high order formula) with stepsize control<a name="line.27"></a>
31 <FONT color="green">028</FONT> * (and automatic step initialization) and continuous output. This<a name="line.28"></a>
32 <FONT color="green">029</FONT> * method uses 12 functions evaluations per step for integration and 4<a name="line.29"></a>
33 <FONT color="green">030</FONT> * evaluations for interpolation. However, since the first<a name="line.30"></a>
34 <FONT color="green">031</FONT> * interpolation evaluation is the same as the first integration<a name="line.31"></a>
35 <FONT color="green">032</FONT> * evaluation of the next step, we have included it in the integrator<a name="line.32"></a>
36 <FONT color="green">033</FONT> * rather than in the interpolator and specified the method was an<a name="line.33"></a>
37 <FONT color="green">034</FONT> * &lt;i&gt;fsal&lt;/i&gt;. Hence, despite we have 13 stages here, the cost is<a name="line.34"></a>
38 <FONT color="green">035</FONT> * really 12 evaluations per step even if no interpolation is done,<a name="line.35"></a>
39 <FONT color="green">036</FONT> * and the overcost of interpolation is only 3 evaluations.&lt;/p&gt;<a name="line.36"></a>
40 <FONT color="green">037</FONT> *<a name="line.37"></a>
41 <FONT color="green">038</FONT> * &lt;p&gt;This method is based on an 8(6) method by Dormand and Prince<a name="line.38"></a>
42 <FONT color="green">039</FONT> * (i.e. order 8 for the integration and order 6 for error estimation)<a name="line.39"></a>
43 <FONT color="green">040</FONT> * modified by Hairer and Wanner to use a 5th order error estimator<a name="line.40"></a>
44 <FONT color="green">041</FONT> * with 3rd order correction. This modification was introduced because<a name="line.41"></a>
45 <FONT color="green">042</FONT> * the original method failed in some cases (wrong steps can be<a name="line.42"></a>
46 <FONT color="green">043</FONT> * accepted when step size is too large, for example in the<a name="line.43"></a>
47 <FONT color="green">044</FONT> * Brusselator problem) and also had &lt;i&gt;severe difficulties when<a name="line.44"></a>
48 <FONT color="green">045</FONT> * applied to problems with discontinuities&lt;/i&gt;. This modification is<a name="line.45"></a>
49 <FONT color="green">046</FONT> * explained in the second edition of the first volume (Nonstiff<a name="line.46"></a>
50 <FONT color="green">047</FONT> * Problems) of the reference book by Hairer, Norsett and Wanner:<a name="line.47"></a>
51 <FONT color="green">048</FONT> * &lt;i&gt;Solving Ordinary Differential Equations&lt;/i&gt; (Springer-Verlag,<a name="line.48"></a>
52 <FONT color="green">049</FONT> * ISBN 3-540-56670-8).&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> * @version $Revision: 810196 $ $Date: 2009-09-01 15:47:46 -0400 (Tue, 01 Sep 2009) $<a name="line.51"></a>
55 <FONT color="green">052</FONT> * @since 1.2<a name="line.52"></a>
56 <FONT color="green">053</FONT> */<a name="line.53"></a>
57 <FONT color="green">054</FONT> <a name="line.54"></a>
58 <FONT color="green">055</FONT> public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {<a name="line.55"></a>
59 <FONT color="green">056</FONT> <a name="line.56"></a>
60 <FONT color="green">057</FONT> /** Integrator method name. */<a name="line.57"></a>
61 <FONT color="green">058</FONT> private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";<a name="line.58"></a>
62 <FONT color="green">059</FONT> <a name="line.59"></a>
63 <FONT color="green">060</FONT> /** Time steps Butcher array. */<a name="line.60"></a>
64 <FONT color="green">061</FONT> private static final double[] STATIC_C = {<a name="line.61"></a>
65 <FONT color="green">062</FONT> (12.0 - 2.0 * Math.sqrt(6.0)) / 135.0, (6.0 - Math.sqrt(6.0)) / 45.0, (6.0 - Math.sqrt(6.0)) / 30.0,<a name="line.62"></a>
66 <FONT color="green">063</FONT> (6.0 + Math.sqrt(6.0)) / 30.0, 1.0/3.0, 1.0/4.0, 4.0/13.0, 127.0/195.0, 3.0/5.0,<a name="line.63"></a>
67 <FONT color="green">064</FONT> 6.0/7.0, 1.0, 1.0<a name="line.64"></a>
68 <FONT color="green">065</FONT> };<a name="line.65"></a>
69 <FONT color="green">066</FONT> <a name="line.66"></a>
70 <FONT color="green">067</FONT> /** Internal weights Butcher array. */<a name="line.67"></a>
71 <FONT color="green">068</FONT> private static final double[][] STATIC_A = {<a name="line.68"></a>
72 <FONT color="green">069</FONT> <a name="line.69"></a>
73 <FONT color="green">070</FONT> // k2<a name="line.70"></a>
74 <FONT color="green">071</FONT> {(12.0 - 2.0 * Math.sqrt(6.0)) / 135.0},<a name="line.71"></a>
75 <FONT color="green">072</FONT> <a name="line.72"></a>
76 <FONT color="green">073</FONT> // k3<a name="line.73"></a>
77 <FONT color="green">074</FONT> {(6.0 - Math.sqrt(6.0)) / 180.0, (6.0 - Math.sqrt(6.0)) / 60.0},<a name="line.74"></a>
78 <FONT color="green">075</FONT> <a name="line.75"></a>
79 <FONT color="green">076</FONT> // k4<a name="line.76"></a>
80 <FONT color="green">077</FONT> {(6.0 - Math.sqrt(6.0)) / 120.0, 0.0, (6.0 - Math.sqrt(6.0)) / 40.0},<a name="line.77"></a>
81 <FONT color="green">078</FONT> <a name="line.78"></a>
82 <FONT color="green">079</FONT> // k5<a name="line.79"></a>
83 <FONT color="green">080</FONT> {(462.0 + 107.0 * Math.sqrt(6.0)) / 3000.0, 0.0,<a name="line.80"></a>
84 <FONT color="green">081</FONT> (-402.0 - 197.0 * Math.sqrt(6.0)) / 1000.0, (168.0 + 73.0 * Math.sqrt(6.0)) / 375.0},<a name="line.81"></a>
85 <FONT color="green">082</FONT> <a name="line.82"></a>
86 <FONT color="green">083</FONT> // k6<a name="line.83"></a>
87 <FONT color="green">084</FONT> {1.0 / 27.0, 0.0, 0.0, (16.0 + Math.sqrt(6.0)) / 108.0, (16.0 - Math.sqrt(6.0)) / 108.0},<a name="line.84"></a>
88 <FONT color="green">085</FONT> <a name="line.85"></a>
89 <FONT color="green">086</FONT> // k7<a name="line.86"></a>
90 <FONT color="green">087</FONT> {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * Math.sqrt(6.0)) / 1024.0,<a name="line.87"></a>
91 <FONT color="green">088</FONT> (118.0 - 23.0 * Math.sqrt(6.0)) / 1024.0, -9.0 / 512.0},<a name="line.88"></a>
92 <FONT color="green">089</FONT> <a name="line.89"></a>
93 <FONT color="green">090</FONT> // k8<a name="line.90"></a>
94 <FONT color="green">091</FONT> {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * Math.sqrt(6.0)) / 371293.0,<a name="line.91"></a>
95 <FONT color="green">092</FONT> (51544.0 - 4784.0 * Math.sqrt(6.0)) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0},<a name="line.92"></a>
96 <FONT color="green">093</FONT> <a name="line.93"></a>
97 <FONT color="green">094</FONT> // k9<a name="line.94"></a>
98 <FONT color="green">095</FONT> {58656157643.0 / 93983540625.0, 0.0, 0.0,<a name="line.95"></a>
99 <FONT color="green">096</FONT> (-1324889724104.0 - 318801444819.0 * Math.sqrt(6.0)) / 626556937500.0,<a name="line.96"></a>
100 <FONT color="green">097</FONT> (-1324889724104.0 + 318801444819.0 * Math.sqrt(6.0)) / 626556937500.0,<a name="line.97"></a>
101 <FONT color="green">098</FONT> 96044563816.0 / 3480871875.0, 5682451879168.0 / 281950621875.0,<a name="line.98"></a>
102 <FONT color="green">099</FONT> -165125654.0 / 3796875.0},<a name="line.99"></a>
103 <FONT color="green">100</FONT> <a name="line.100"></a>
104 <FONT color="green">101</FONT> // k10<a name="line.101"></a>
105 <FONT color="green">102</FONT> {8909899.0 / 18653125.0, 0.0, 0.0,<a name="line.102"></a>
106 <FONT color="green">103</FONT> (-4521408.0 - 1137963.0 * Math.sqrt(6.0)) / 2937500.0,<a name="line.103"></a>
107 <FONT color="green">104</FONT> (-4521408.0 + 1137963.0 * Math.sqrt(6.0)) / 2937500.0,<a name="line.104"></a>
108 <FONT color="green">105</FONT> 96663078.0 / 4553125.0, 2107245056.0 / 137915625.0,<a name="line.105"></a>
109 <FONT color="green">106</FONT> -4913652016.0 / 147609375.0, -78894270.0 / 3880452869.0},<a name="line.106"></a>
110 <FONT color="green">107</FONT> <a name="line.107"></a>
111 <FONT color="green">108</FONT> // k11<a name="line.108"></a>
112 <FONT color="green">109</FONT> {-20401265806.0 / 21769653311.0, 0.0, 0.0,<a name="line.109"></a>
113 <FONT color="green">110</FONT> (354216.0 + 94326.0 * Math.sqrt(6.0)) / 112847.0,<a name="line.110"></a>
114 <FONT color="green">111</FONT> (354216.0 - 94326.0 * Math.sqrt(6.0)) / 112847.0,<a name="line.111"></a>
115 <FONT color="green">112</FONT> -43306765128.0 / 5313852383.0, -20866708358144.0 / 1126708119789.0,<a name="line.112"></a>
116 <FONT color="green">113</FONT> 14886003438020.0 / 654632330667.0, 35290686222309375.0 / 14152473387134411.0,<a name="line.113"></a>
117 <FONT color="green">114</FONT> -1477884375.0 / 485066827.0},<a name="line.114"></a>
118 <FONT color="green">115</FONT> <a name="line.115"></a>
119 <FONT color="green">116</FONT> // k12<a name="line.116"></a>
120 <FONT color="green">117</FONT> {39815761.0 / 17514443.0, 0.0, 0.0,<a name="line.117"></a>
121 <FONT color="green">118</FONT> (-3457480.0 - 960905.0 * Math.sqrt(6.0)) / 551636.0,<a name="line.118"></a>
122 <FONT color="green">119</FONT> (-3457480.0 + 960905.0 * Math.sqrt(6.0)) / 551636.0,<a name="line.119"></a>
123 <FONT color="green">120</FONT> -844554132.0 / 47026969.0, 8444996352.0 / 302158619.0,<a name="line.120"></a>
124 <FONT color="green">121</FONT> -2509602342.0 / 877790785.0, -28388795297996250.0 / 3199510091356783.0,<a name="line.121"></a>
125 <FONT color="green">122</FONT> 226716250.0 / 18341897.0, 1371316744.0 / 2131383595.0},<a name="line.122"></a>
126 <FONT color="green">123</FONT> <a name="line.123"></a>
127 <FONT color="green">124</FONT> // k13 should be for interpolation only, but since it is the same<a name="line.124"></a>
128 <FONT color="green">125</FONT> // stage as the first evaluation of the next step, we perform it<a name="line.125"></a>
129 <FONT color="green">126</FONT> // here at no cost by specifying this is an fsal method<a name="line.126"></a>
130 <FONT color="green">127</FONT> {104257.0/1920240.0, 0.0, 0.0, 0.0, 0.0, 3399327.0/763840.0,<a name="line.127"></a>
131 <FONT color="green">128</FONT> 66578432.0/35198415.0, -1674902723.0/288716400.0,<a name="line.128"></a>
132 <FONT color="green">129</FONT> 54980371265625.0/176692375811392.0, -734375.0/4826304.0,<a name="line.129"></a>
133 <FONT color="green">130</FONT> 171414593.0/851261400.0, 137909.0/3084480.0}<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> /** Propagation weights Butcher array. */<a name="line.134"></a>
138 <FONT color="green">135</FONT> private static final double[] STATIC_B = {<a name="line.135"></a>
139 <FONT color="green">136</FONT> 104257.0/1920240.0,<a name="line.136"></a>
140 <FONT color="green">137</FONT> 0.0,<a name="line.137"></a>
141 <FONT color="green">138</FONT> 0.0,<a name="line.138"></a>
142 <FONT color="green">139</FONT> 0.0,<a name="line.139"></a>
143 <FONT color="green">140</FONT> 0.0,<a name="line.140"></a>
144 <FONT color="green">141</FONT> 3399327.0/763840.0,<a name="line.141"></a>
145 <FONT color="green">142</FONT> 66578432.0/35198415.0,<a name="line.142"></a>
146 <FONT color="green">143</FONT> -1674902723.0/288716400.0,<a name="line.143"></a>
147 <FONT color="green">144</FONT> 54980371265625.0/176692375811392.0,<a name="line.144"></a>
148 <FONT color="green">145</FONT> -734375.0/4826304.0,<a name="line.145"></a>
149 <FONT color="green">146</FONT> 171414593.0/851261400.0,<a name="line.146"></a>
150 <FONT color="green">147</FONT> 137909.0/3084480.0,<a name="line.147"></a>
151 <FONT color="green">148</FONT> 0.0<a name="line.148"></a>
152 <FONT color="green">149</FONT> };<a name="line.149"></a>
153 <FONT color="green">150</FONT> <a name="line.150"></a>
154 <FONT color="green">151</FONT> /** First error weights array, element 1. */<a name="line.151"></a>
155 <FONT color="green">152</FONT> private static final double E1_01 = 116092271.0 / 8848465920.0;<a name="line.152"></a>
156 <FONT color="green">153</FONT> <a name="line.153"></a>
157 <FONT color="green">154</FONT> // elements 2 to 5 are zero, so they are neither stored nor used<a name="line.154"></a>
158 <FONT color="green">155</FONT> <a name="line.155"></a>
159 <FONT color="green">156</FONT> /** First error weights array, element 6. */<a name="line.156"></a>
160 <FONT color="green">157</FONT> private static final double E1_06 = -1871647.0 / 1527680.0;<a name="line.157"></a>
161 <FONT color="green">158</FONT> <a name="line.158"></a>
162 <FONT color="green">159</FONT> /** First error weights array, element 7. */<a name="line.159"></a>
163 <FONT color="green">160</FONT> private static final double E1_07 = -69799717.0 / 140793660.0;<a name="line.160"></a>
164 <FONT color="green">161</FONT> <a name="line.161"></a>
165 <FONT color="green">162</FONT> /** First error weights array, element 8. */<a name="line.162"></a>
166 <FONT color="green">163</FONT> private static final double E1_08 = 1230164450203.0 / 739113984000.0;<a name="line.163"></a>
167 <FONT color="green">164</FONT> <a name="line.164"></a>
168 <FONT color="green">165</FONT> /** First error weights array, element 9. */<a name="line.165"></a>
169 <FONT color="green">166</FONT> private static final double E1_09 = -1980813971228885.0 / 5654156025964544.0;<a name="line.166"></a>
170 <FONT color="green">167</FONT> <a name="line.167"></a>
171 <FONT color="green">168</FONT> /** First error weights array, element 10. */<a name="line.168"></a>
172 <FONT color="green">169</FONT> private static final double E1_10 = 464500805.0 / 1389975552.0;<a name="line.169"></a>
173 <FONT color="green">170</FONT> <a name="line.170"></a>
174 <FONT color="green">171</FONT> /** First error weights array, element 11. */<a name="line.171"></a>
175 <FONT color="green">172</FONT> private static final double E1_11 = 1606764981773.0 / 19613062656000.0;<a name="line.172"></a>
176 <FONT color="green">173</FONT> <a name="line.173"></a>
177 <FONT color="green">174</FONT> /** First error weights array, element 12. */<a name="line.174"></a>
178 <FONT color="green">175</FONT> private static final double E1_12 = -137909.0 / 6168960.0;<a name="line.175"></a>
179 <FONT color="green">176</FONT> <a name="line.176"></a>
180 <FONT color="green">177</FONT> <a name="line.177"></a>
181 <FONT color="green">178</FONT> /** Second error weights array, element 1. */<a name="line.178"></a>
182 <FONT color="green">179</FONT> private static final double E2_01 = -364463.0 / 1920240.0;<a name="line.179"></a>
183 <FONT color="green">180</FONT> <a name="line.180"></a>
184 <FONT color="green">181</FONT> // elements 2 to 5 are zero, so they are neither stored nor used<a name="line.181"></a>
185 <FONT color="green">182</FONT> <a name="line.182"></a>
186 <FONT color="green">183</FONT> /** Second error weights array, element 6. */<a name="line.183"></a>
187 <FONT color="green">184</FONT> private static final double E2_06 = 3399327.0 / 763840.0;<a name="line.184"></a>
188 <FONT color="green">185</FONT> <a name="line.185"></a>
189 <FONT color="green">186</FONT> /** Second error weights array, element 7. */<a name="line.186"></a>
190 <FONT color="green">187</FONT> private static final double E2_07 = 66578432.0 / 35198415.0;<a name="line.187"></a>
191 <FONT color="green">188</FONT> <a name="line.188"></a>
192 <FONT color="green">189</FONT> /** Second error weights array, element 8. */<a name="line.189"></a>
193 <FONT color="green">190</FONT> private static final double E2_08 = -1674902723.0 / 288716400.0;<a name="line.190"></a>
194 <FONT color="green">191</FONT> <a name="line.191"></a>
195 <FONT color="green">192</FONT> /** Second error weights array, element 9. */<a name="line.192"></a>
196 <FONT color="green">193</FONT> private static final double E2_09 = -74684743568175.0 / 176692375811392.0;<a name="line.193"></a>
197 <FONT color="green">194</FONT> <a name="line.194"></a>
198 <FONT color="green">195</FONT> /** Second error weights array, element 10. */<a name="line.195"></a>
199 <FONT color="green">196</FONT> private static final double E2_10 = -734375.0 / 4826304.0;<a name="line.196"></a>
200 <FONT color="green">197</FONT> <a name="line.197"></a>
201 <FONT color="green">198</FONT> /** Second error weights array, element 11. */<a name="line.198"></a>
202 <FONT color="green">199</FONT> private static final double E2_11 = 171414593.0 / 851261400.0;<a name="line.199"></a>
203 <FONT color="green">200</FONT> <a name="line.200"></a>
204 <FONT color="green">201</FONT> /** Second error weights array, element 12. */<a name="line.201"></a>
205 <FONT color="green">202</FONT> private static final double E2_12 = 69869.0 / 3084480.0;<a name="line.202"></a>
206 <FONT color="green">203</FONT> <a name="line.203"></a>
207 <FONT color="green">204</FONT> /** Simple constructor.<a name="line.204"></a>
208 <FONT color="green">205</FONT> * Build an eighth order Dormand-Prince integrator with the given step bounds<a name="line.205"></a>
209 <FONT color="green">206</FONT> * @param minStep minimal step (must be positive even for backward<a name="line.206"></a>
210 <FONT color="green">207</FONT> * integration), the last step can be smaller than this<a name="line.207"></a>
211 <FONT color="green">208</FONT> * @param maxStep maximal step (must be positive even for backward<a name="line.208"></a>
212 <FONT color="green">209</FONT> * integration)<a name="line.209"></a>
213 <FONT color="green">210</FONT> * @param scalAbsoluteTolerance allowed absolute error<a name="line.210"></a>
214 <FONT color="green">211</FONT> * @param scalRelativeTolerance allowed relative error<a name="line.211"></a>
215 <FONT color="green">212</FONT> */<a name="line.212"></a>
216 <FONT color="green">213</FONT> public DormandPrince853Integrator(final double minStep, final double maxStep,<a name="line.213"></a>
217 <FONT color="green">214</FONT> final double scalAbsoluteTolerance,<a name="line.214"></a>
218 <FONT color="green">215</FONT> final double scalRelativeTolerance) {<a name="line.215"></a>
219 <FONT color="green">216</FONT> super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,<a name="line.216"></a>
220 <FONT color="green">217</FONT> new DormandPrince853StepInterpolator(),<a name="line.217"></a>
221 <FONT color="green">218</FONT> minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);<a name="line.218"></a>
222 <FONT color="green">219</FONT> }<a name="line.219"></a>
223 <FONT color="green">220</FONT> <a name="line.220"></a>
224 <FONT color="green">221</FONT> /** Simple constructor.<a name="line.221"></a>
225 <FONT color="green">222</FONT> * Build an eighth order Dormand-Prince integrator with the given step bounds<a name="line.222"></a>
226 <FONT color="green">223</FONT> * @param minStep minimal step (must be positive even for backward<a name="line.223"></a>
227 <FONT color="green">224</FONT> * integration), the last step can be smaller than this<a name="line.224"></a>
228 <FONT color="green">225</FONT> * @param maxStep maximal step (must be positive even for backward<a name="line.225"></a>
229 <FONT color="green">226</FONT> * integration)<a name="line.226"></a>
230 <FONT color="green">227</FONT> * @param vecAbsoluteTolerance allowed absolute error<a name="line.227"></a>
231 <FONT color="green">228</FONT> * @param vecRelativeTolerance allowed relative error<a name="line.228"></a>
232 <FONT color="green">229</FONT> */<a name="line.229"></a>
233 <FONT color="green">230</FONT> public DormandPrince853Integrator(final double minStep, final double maxStep,<a name="line.230"></a>
234 <FONT color="green">231</FONT> final double[] vecAbsoluteTolerance,<a name="line.231"></a>
235 <FONT color="green">232</FONT> final double[] vecRelativeTolerance) {<a name="line.232"></a>
236 <FONT color="green">233</FONT> super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,<a name="line.233"></a>
237 <FONT color="green">234</FONT> new DormandPrince853StepInterpolator(),<a name="line.234"></a>
238 <FONT color="green">235</FONT> minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);<a name="line.235"></a>
239 <FONT color="green">236</FONT> }<a name="line.236"></a>
240 <FONT color="green">237</FONT> <a name="line.237"></a>
241 <FONT color="green">238</FONT> /** {@inheritDoc} */<a name="line.238"></a>
242 <FONT color="green">239</FONT> @Override<a name="line.239"></a>
243 <FONT color="green">240</FONT> public int getOrder() {<a name="line.240"></a>
244 <FONT color="green">241</FONT> return 8;<a name="line.241"></a>
245 <FONT color="green">242</FONT> }<a name="line.242"></a>
246 <FONT color="green">243</FONT> <a name="line.243"></a>
247 <FONT color="green">244</FONT> /** {@inheritDoc} */<a name="line.244"></a>
248 <FONT color="green">245</FONT> @Override<a name="line.245"></a>
249 <FONT color="green">246</FONT> protected double estimateError(final double[][] yDotK,<a name="line.246"></a>
250 <FONT color="green">247</FONT> final double[] y0, final double[] y1,<a name="line.247"></a>
251 <FONT color="green">248</FONT> final double h) {<a name="line.248"></a>
252 <FONT color="green">249</FONT> double error1 = 0;<a name="line.249"></a>
253 <FONT color="green">250</FONT> double error2 = 0;<a name="line.250"></a>
254 <FONT color="green">251</FONT> <a name="line.251"></a>
255 <FONT color="green">252</FONT> for (int j = 0; j &lt; y0.length; ++j) {<a name="line.252"></a>
256 <FONT color="green">253</FONT> final double errSum1 = E1_01 * yDotK[0][j] + E1_06 * yDotK[5][j] +<a name="line.253"></a>
257 <FONT color="green">254</FONT> E1_07 * yDotK[6][j] + E1_08 * yDotK[7][j] +<a name="line.254"></a>
258 <FONT color="green">255</FONT> E1_09 * yDotK[8][j] + E1_10 * yDotK[9][j] +<a name="line.255"></a>
259 <FONT color="green">256</FONT> E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j];<a name="line.256"></a>
260 <FONT color="green">257</FONT> final double errSum2 = E2_01 * yDotK[0][j] + E2_06 * yDotK[5][j] +<a name="line.257"></a>
261 <FONT color="green">258</FONT> E2_07 * yDotK[6][j] + E2_08 * yDotK[7][j] +<a name="line.258"></a>
262 <FONT color="green">259</FONT> E2_09 * yDotK[8][j] + E2_10 * yDotK[9][j] +<a name="line.259"></a>
263 <FONT color="green">260</FONT> E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j];<a name="line.260"></a>
264 <FONT color="green">261</FONT> <a name="line.261"></a>
265 <FONT color="green">262</FONT> final double yScale = Math.max(Math.abs(y0[j]), Math.abs(y1[j]));<a name="line.262"></a>
266 <FONT color="green">263</FONT> final double tol = (vecAbsoluteTolerance == null) ?<a name="line.263"></a>
267 <FONT color="green">264</FONT> (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :<a name="line.264"></a>
268 <FONT color="green">265</FONT> (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);<a name="line.265"></a>
269 <FONT color="green">266</FONT> final double ratio1 = errSum1 / tol;<a name="line.266"></a>
270 <FONT color="green">267</FONT> error1 += ratio1 * ratio1;<a name="line.267"></a>
271 <FONT color="green">268</FONT> final double ratio2 = errSum2 / tol;<a name="line.268"></a>
272 <FONT color="green">269</FONT> error2 += ratio2 * ratio2;<a name="line.269"></a>
273 <FONT color="green">270</FONT> }<a name="line.270"></a>
274 <FONT color="green">271</FONT> <a name="line.271"></a>
275 <FONT color="green">272</FONT> double den = error1 + 0.01 * error2;<a name="line.272"></a>
276 <FONT color="green">273</FONT> if (den &lt;= 0.0) {<a name="line.273"></a>
277 <FONT color="green">274</FONT> den = 1.0;<a name="line.274"></a>
278 <FONT color="green">275</FONT> }<a name="line.275"></a>
279 <FONT color="green">276</FONT> <a name="line.276"></a>
280 <FONT color="green">277</FONT> return Math.abs(h) * error1 / Math.sqrt(y0.length * den);<a name="line.277"></a>
281 <FONT color="green">278</FONT> <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> }<a name="line.281"></a>
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
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 </PRE>
346 </BODY>
347 </HTML>