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

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/analysis/solvers/BrentSolver.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> package org.apache.commons.math.analysis.solvers;<a name="line.17"></a>
21 <FONT color="green">018</FONT> <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.FunctionEvaluationException;<a name="line.20"></a>
24 <FONT color="green">021</FONT> import org.apache.commons.math.MathRuntimeException;<a name="line.21"></a>
25 <FONT color="green">022</FONT> import org.apache.commons.math.MaxIterationsExceededException;<a name="line.22"></a>
26 <FONT color="green">023</FONT> import org.apache.commons.math.analysis.UnivariateRealFunction;<a name="line.23"></a>
27 <FONT color="green">024</FONT> <a name="line.24"></a>
28 <FONT color="green">025</FONT> /**<a name="line.25"></a>
29 <FONT color="green">026</FONT> * Implements the &lt;a href="http://mathworld.wolfram.com/BrentsMethod.html"&gt;<a name="line.26"></a>
30 <FONT color="green">027</FONT> * Brent algorithm&lt;/a&gt; for finding zeros of real univariate functions.<a name="line.27"></a>
31 <FONT color="green">028</FONT> * &lt;p&gt;<a name="line.28"></a>
32 <FONT color="green">029</FONT> * The function should be continuous but not necessarily smooth.&lt;/p&gt;<a name="line.29"></a>
33 <FONT color="green">030</FONT> *<a name="line.30"></a>
34 <FONT color="green">031</FONT> * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $<a name="line.31"></a>
35 <FONT color="green">032</FONT> */<a name="line.32"></a>
36 <FONT color="green">033</FONT> public class BrentSolver extends UnivariateRealSolverImpl {<a name="line.33"></a>
37 <FONT color="green">034</FONT> <a name="line.34"></a>
38 <FONT color="green">035</FONT> /**<a name="line.35"></a>
39 <FONT color="green">036</FONT> * Default absolute accuracy<a name="line.36"></a>
40 <FONT color="green">037</FONT> * @since 2.1<a name="line.37"></a>
41 <FONT color="green">038</FONT> */<a name="line.38"></a>
42 <FONT color="green">039</FONT> public static final double DEFAULT_ABSOLUTE_ACCURACY = 1E-6;<a name="line.39"></a>
43 <FONT color="green">040</FONT> <a name="line.40"></a>
44 <FONT color="green">041</FONT> /** Default maximum number of iterations<a name="line.41"></a>
45 <FONT color="green">042</FONT> * @since 2.1<a name="line.42"></a>
46 <FONT color="green">043</FONT> */<a name="line.43"></a>
47 <FONT color="green">044</FONT> public static final int DEFAULT_MAXIMUM_ITERATIONS = 100;<a name="line.44"></a>
48 <FONT color="green">045</FONT> <a name="line.45"></a>
49 <FONT color="green">046</FONT> /** Error message for non-bracketing interval. */<a name="line.46"></a>
50 <FONT color="green">047</FONT> private static final String NON_BRACKETING_MESSAGE =<a name="line.47"></a>
51 <FONT color="green">048</FONT> "function values at endpoints do not have different signs. " +<a name="line.48"></a>
52 <FONT color="green">049</FONT> "Endpoints: [{0}, {1}], Values: [{2}, {3}]";<a name="line.49"></a>
53 <FONT color="green">050</FONT> <a name="line.50"></a>
54 <FONT color="green">051</FONT> /** Serializable version identifier */<a name="line.51"></a>
55 <FONT color="green">052</FONT> private static final long serialVersionUID = 7694577816772532779L;<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> * Construct a solver for the given function.<a name="line.55"></a>
59 <FONT color="green">056</FONT> *<a name="line.56"></a>
60 <FONT color="green">057</FONT> * @param f function to solve.<a name="line.57"></a>
61 <FONT color="green">058</FONT> * @deprecated as of 2.0 the function to solve is passed as an argument<a name="line.58"></a>
62 <FONT color="green">059</FONT> * to the {@link #solve(UnivariateRealFunction, double, double)} or<a name="line.59"></a>
63 <FONT color="green">060</FONT> * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}<a name="line.60"></a>
64 <FONT color="green">061</FONT> * method.<a name="line.61"></a>
65 <FONT color="green">062</FONT> */<a name="line.62"></a>
66 <FONT color="green">063</FONT> @Deprecated<a name="line.63"></a>
67 <FONT color="green">064</FONT> public BrentSolver(UnivariateRealFunction f) {<a name="line.64"></a>
68 <FONT color="green">065</FONT> super(f, DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);<a name="line.65"></a>
69 <FONT color="green">066</FONT> }<a name="line.66"></a>
70 <FONT color="green">067</FONT> <a name="line.67"></a>
71 <FONT color="green">068</FONT> /**<a name="line.68"></a>
72 <FONT color="green">069</FONT> * Construct a solver with default properties.<a name="line.69"></a>
73 <FONT color="green">070</FONT> */<a name="line.70"></a>
74 <FONT color="green">071</FONT> public BrentSolver() {<a name="line.71"></a>
75 <FONT color="green">072</FONT> super(DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);<a name="line.72"></a>
76 <FONT color="green">073</FONT> }<a name="line.73"></a>
77 <FONT color="green">074</FONT> <a name="line.74"></a>
78 <FONT color="green">075</FONT> /**<a name="line.75"></a>
79 <FONT color="green">076</FONT> * Construct a solver with the given absolute accuracy.<a name="line.76"></a>
80 <FONT color="green">077</FONT> *<a name="line.77"></a>
81 <FONT color="green">078</FONT> * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver<a name="line.78"></a>
82 <FONT color="green">079</FONT> * @since 2.1<a name="line.79"></a>
83 <FONT color="green">080</FONT> */<a name="line.80"></a>
84 <FONT color="green">081</FONT> public BrentSolver(double absoluteAccuracy) {<a name="line.81"></a>
85 <FONT color="green">082</FONT> super(DEFAULT_MAXIMUM_ITERATIONS, absoluteAccuracy);<a name="line.82"></a>
86 <FONT color="green">083</FONT> }<a name="line.83"></a>
87 <FONT color="green">084</FONT> <a name="line.84"></a>
88 <FONT color="green">085</FONT> /**<a name="line.85"></a>
89 <FONT color="green">086</FONT> * Contstruct a solver with the given maximum iterations and absolute accuracy.<a name="line.86"></a>
90 <FONT color="green">087</FONT> *<a name="line.87"></a>
91 <FONT color="green">088</FONT> * @param maximumIterations maximum number of iterations<a name="line.88"></a>
92 <FONT color="green">089</FONT> * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver<a name="line.89"></a>
93 <FONT color="green">090</FONT> * @since 2.1<a name="line.90"></a>
94 <FONT color="green">091</FONT> */<a name="line.91"></a>
95 <FONT color="green">092</FONT> public BrentSolver(int maximumIterations, double absoluteAccuracy) {<a name="line.92"></a>
96 <FONT color="green">093</FONT> super(maximumIterations, absoluteAccuracy);<a name="line.93"></a>
97 <FONT color="green">094</FONT> }<a name="line.94"></a>
98 <FONT color="green">095</FONT> <a name="line.95"></a>
99 <FONT color="green">096</FONT> /** {@inheritDoc} */<a name="line.96"></a>
100 <FONT color="green">097</FONT> @Deprecated<a name="line.97"></a>
101 <FONT color="green">098</FONT> public double solve(double min, double max)<a name="line.98"></a>
102 <FONT color="green">099</FONT> throws MaxIterationsExceededException, FunctionEvaluationException {<a name="line.99"></a>
103 <FONT color="green">100</FONT> return solve(f, min, max);<a name="line.100"></a>
104 <FONT color="green">101</FONT> }<a name="line.101"></a>
105 <FONT color="green">102</FONT> <a name="line.102"></a>
106 <FONT color="green">103</FONT> /** {@inheritDoc} */<a name="line.103"></a>
107 <FONT color="green">104</FONT> @Deprecated<a name="line.104"></a>
108 <FONT color="green">105</FONT> public double solve(double min, double max, double initial)<a name="line.105"></a>
109 <FONT color="green">106</FONT> throws MaxIterationsExceededException, FunctionEvaluationException {<a name="line.106"></a>
110 <FONT color="green">107</FONT> return solve(f, min, max, initial);<a name="line.107"></a>
111 <FONT color="green">108</FONT> }<a name="line.108"></a>
112 <FONT color="green">109</FONT> <a name="line.109"></a>
113 <FONT color="green">110</FONT> /**<a name="line.110"></a>
114 <FONT color="green">111</FONT> * Find a zero in the given interval with an initial guess.<a name="line.111"></a>
115 <FONT color="green">112</FONT> * &lt;p&gt;Throws &lt;code&gt;IllegalArgumentException&lt;/code&gt; if the values of the<a name="line.112"></a>
116 <FONT color="green">113</FONT> * function at the three points have the same sign (note that it is<a name="line.113"></a>
117 <FONT color="green">114</FONT> * allowed to have endpoints with the same sign if the initial point has<a name="line.114"></a>
118 <FONT color="green">115</FONT> * opposite sign function-wise).&lt;/p&gt;<a name="line.115"></a>
119 <FONT color="green">116</FONT> *<a name="line.116"></a>
120 <FONT color="green">117</FONT> * @param f function to solve.<a name="line.117"></a>
121 <FONT color="green">118</FONT> * @param min the lower bound for the interval.<a name="line.118"></a>
122 <FONT color="green">119</FONT> * @param max the upper bound for the interval.<a name="line.119"></a>
123 <FONT color="green">120</FONT> * @param initial the start value to use (must be set to min if no<a name="line.120"></a>
124 <FONT color="green">121</FONT> * initial point is known).<a name="line.121"></a>
125 <FONT color="green">122</FONT> * @return the value where the function is zero<a name="line.122"></a>
126 <FONT color="green">123</FONT> * @throws MaxIterationsExceededException the maximum iteration count<a name="line.123"></a>
127 <FONT color="green">124</FONT> * is exceeded<a name="line.124"></a>
128 <FONT color="green">125</FONT> * @throws FunctionEvaluationException if an error occurs evaluating<a name="line.125"></a>
129 <FONT color="green">126</FONT> * the function<a name="line.126"></a>
130 <FONT color="green">127</FONT> * @throws IllegalArgumentException if initial is not between min and max<a name="line.127"></a>
131 <FONT color="green">128</FONT> * (even if it &lt;em&gt;is&lt;/em&gt; a root)<a name="line.128"></a>
132 <FONT color="green">129</FONT> */<a name="line.129"></a>
133 <FONT color="green">130</FONT> public double solve(final UnivariateRealFunction f,<a name="line.130"></a>
134 <FONT color="green">131</FONT> final double min, final double max, final double initial)<a name="line.131"></a>
135 <FONT color="green">132</FONT> throws MaxIterationsExceededException, FunctionEvaluationException {<a name="line.132"></a>
136 <FONT color="green">133</FONT> <a name="line.133"></a>
137 <FONT color="green">134</FONT> clearResult();<a name="line.134"></a>
138 <FONT color="green">135</FONT> if ((initial &lt; min) || (initial &gt; max)) {<a name="line.135"></a>
139 <FONT color="green">136</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.136"></a>
140 <FONT color="green">137</FONT> "invalid interval, initial value parameters: lower={0}, initial={1}, upper={2}",<a name="line.137"></a>
141 <FONT color="green">138</FONT> min, initial, max);<a name="line.138"></a>
142 <FONT color="green">139</FONT> }<a name="line.139"></a>
143 <FONT color="green">140</FONT> <a name="line.140"></a>
144 <FONT color="green">141</FONT> // return the initial guess if it is good enough<a name="line.141"></a>
145 <FONT color="green">142</FONT> double yInitial = f.value(initial);<a name="line.142"></a>
146 <FONT color="green">143</FONT> if (Math.abs(yInitial) &lt;= functionValueAccuracy) {<a name="line.143"></a>
147 <FONT color="green">144</FONT> setResult(initial, 0);<a name="line.144"></a>
148 <FONT color="green">145</FONT> return result;<a name="line.145"></a>
149 <FONT color="green">146</FONT> }<a name="line.146"></a>
150 <FONT color="green">147</FONT> <a name="line.147"></a>
151 <FONT color="green">148</FONT> // return the first endpoint if it is good enough<a name="line.148"></a>
152 <FONT color="green">149</FONT> double yMin = f.value(min);<a name="line.149"></a>
153 <FONT color="green">150</FONT> if (Math.abs(yMin) &lt;= functionValueAccuracy) {<a name="line.150"></a>
154 <FONT color="green">151</FONT> setResult(min, 0);<a name="line.151"></a>
155 <FONT color="green">152</FONT> return result;<a name="line.152"></a>
156 <FONT color="green">153</FONT> }<a name="line.153"></a>
157 <FONT color="green">154</FONT> <a name="line.154"></a>
158 <FONT color="green">155</FONT> // reduce interval if min and initial bracket the root<a name="line.155"></a>
159 <FONT color="green">156</FONT> if (yInitial * yMin &lt; 0) {<a name="line.156"></a>
160 <FONT color="green">157</FONT> return solve(f, min, yMin, initial, yInitial, min, yMin);<a name="line.157"></a>
161 <FONT color="green">158</FONT> }<a name="line.158"></a>
162 <FONT color="green">159</FONT> <a name="line.159"></a>
163 <FONT color="green">160</FONT> // return the second endpoint if it is good enough<a name="line.160"></a>
164 <FONT color="green">161</FONT> double yMax = f.value(max);<a name="line.161"></a>
165 <FONT color="green">162</FONT> if (Math.abs(yMax) &lt;= functionValueAccuracy) {<a name="line.162"></a>
166 <FONT color="green">163</FONT> setResult(max, 0);<a name="line.163"></a>
167 <FONT color="green">164</FONT> return result;<a name="line.164"></a>
168 <FONT color="green">165</FONT> }<a name="line.165"></a>
169 <FONT color="green">166</FONT> <a name="line.166"></a>
170 <FONT color="green">167</FONT> // reduce interval if initial and max bracket the root<a name="line.167"></a>
171 <FONT color="green">168</FONT> if (yInitial * yMax &lt; 0) {<a name="line.168"></a>
172 <FONT color="green">169</FONT> return solve(f, initial, yInitial, max, yMax, initial, yInitial);<a name="line.169"></a>
173 <FONT color="green">170</FONT> }<a name="line.170"></a>
174 <FONT color="green">171</FONT> <a name="line.171"></a>
175 <FONT color="green">172</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.172"></a>
176 <FONT color="green">173</FONT> NON_BRACKETING_MESSAGE, min, max, yMin, yMax);<a name="line.173"></a>
177 <FONT color="green">174</FONT> <a name="line.174"></a>
178 <FONT color="green">175</FONT> }<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> * Find a zero in the given interval.<a name="line.178"></a>
182 <FONT color="green">179</FONT> * &lt;p&gt;<a name="line.179"></a>
183 <FONT color="green">180</FONT> * Requires that the values of the function at the endpoints have opposite<a name="line.180"></a>
184 <FONT color="green">181</FONT> * signs. An &lt;code&gt;IllegalArgumentException&lt;/code&gt; is thrown if this is not<a name="line.181"></a>
185 <FONT color="green">182</FONT> * the case.&lt;/p&gt;<a name="line.182"></a>
186 <FONT color="green">183</FONT> *<a name="line.183"></a>
187 <FONT color="green">184</FONT> * @param f the function to solve<a name="line.184"></a>
188 <FONT color="green">185</FONT> * @param min the lower bound for the interval.<a name="line.185"></a>
189 <FONT color="green">186</FONT> * @param max the upper bound for the interval.<a name="line.186"></a>
190 <FONT color="green">187</FONT> * @return the value where the function is zero<a name="line.187"></a>
191 <FONT color="green">188</FONT> * @throws MaxIterationsExceededException if the maximum iteration count is exceeded<a name="line.188"></a>
192 <FONT color="green">189</FONT> * @throws FunctionEvaluationException if an error occurs evaluating the<a name="line.189"></a>
193 <FONT color="green">190</FONT> * function<a name="line.190"></a>
194 <FONT color="green">191</FONT> * @throws IllegalArgumentException if min is not less than max or the<a name="line.191"></a>
195 <FONT color="green">192</FONT> * signs of the values of the function at the endpoints are not opposites<a name="line.192"></a>
196 <FONT color="green">193</FONT> */<a name="line.193"></a>
197 <FONT color="green">194</FONT> public double solve(final UnivariateRealFunction f,<a name="line.194"></a>
198 <FONT color="green">195</FONT> final double min, final double max)<a name="line.195"></a>
199 <FONT color="green">196</FONT> throws MaxIterationsExceededException,<a name="line.196"></a>
200 <FONT color="green">197</FONT> FunctionEvaluationException {<a name="line.197"></a>
201 <FONT color="green">198</FONT> <a name="line.198"></a>
202 <FONT color="green">199</FONT> clearResult();<a name="line.199"></a>
203 <FONT color="green">200</FONT> verifyInterval(min, max);<a name="line.200"></a>
204 <FONT color="green">201</FONT> <a name="line.201"></a>
205 <FONT color="green">202</FONT> double ret = Double.NaN;<a name="line.202"></a>
206 <FONT color="green">203</FONT> <a name="line.203"></a>
207 <FONT color="green">204</FONT> double yMin = f.value(min);<a name="line.204"></a>
208 <FONT color="green">205</FONT> double yMax = f.value(max);<a name="line.205"></a>
209 <FONT color="green">206</FONT> <a name="line.206"></a>
210 <FONT color="green">207</FONT> // Verify bracketing<a name="line.207"></a>
211 <FONT color="green">208</FONT> double sign = yMin * yMax;<a name="line.208"></a>
212 <FONT color="green">209</FONT> if (sign &gt; 0) {<a name="line.209"></a>
213 <FONT color="green">210</FONT> // check if either value is close to a zero<a name="line.210"></a>
214 <FONT color="green">211</FONT> if (Math.abs(yMin) &lt;= functionValueAccuracy) {<a name="line.211"></a>
215 <FONT color="green">212</FONT> setResult(min, 0);<a name="line.212"></a>
216 <FONT color="green">213</FONT> ret = min;<a name="line.213"></a>
217 <FONT color="green">214</FONT> } else if (Math.abs(yMax) &lt;= functionValueAccuracy) {<a name="line.214"></a>
218 <FONT color="green">215</FONT> setResult(max, 0);<a name="line.215"></a>
219 <FONT color="green">216</FONT> ret = max;<a name="line.216"></a>
220 <FONT color="green">217</FONT> } else {<a name="line.217"></a>
221 <FONT color="green">218</FONT> // neither value is close to zero and min and max do not bracket root.<a name="line.218"></a>
222 <FONT color="green">219</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.219"></a>
223 <FONT color="green">220</FONT> NON_BRACKETING_MESSAGE, min, max, yMin, yMax);<a name="line.220"></a>
224 <FONT color="green">221</FONT> }<a name="line.221"></a>
225 <FONT color="green">222</FONT> } else if (sign &lt; 0){<a name="line.222"></a>
226 <FONT color="green">223</FONT> // solve using only the first endpoint as initial guess<a name="line.223"></a>
227 <FONT color="green">224</FONT> ret = solve(f, min, yMin, max, yMax, min, yMin);<a name="line.224"></a>
228 <FONT color="green">225</FONT> } else {<a name="line.225"></a>
229 <FONT color="green">226</FONT> // either min or max is a root<a name="line.226"></a>
230 <FONT color="green">227</FONT> if (yMin == 0.0) {<a name="line.227"></a>
231 <FONT color="green">228</FONT> ret = min;<a name="line.228"></a>
232 <FONT color="green">229</FONT> } else {<a name="line.229"></a>
233 <FONT color="green">230</FONT> ret = max;<a name="line.230"></a>
234 <FONT color="green">231</FONT> }<a name="line.231"></a>
235 <FONT color="green">232</FONT> }<a name="line.232"></a>
236 <FONT color="green">233</FONT> <a name="line.233"></a>
237 <FONT color="green">234</FONT> return ret;<a name="line.234"></a>
238 <FONT color="green">235</FONT> }<a name="line.235"></a>
239 <FONT color="green">236</FONT> <a name="line.236"></a>
240 <FONT color="green">237</FONT> /**<a name="line.237"></a>
241 <FONT color="green">238</FONT> * Find a zero starting search according to the three provided points.<a name="line.238"></a>
242 <FONT color="green">239</FONT> * @param f the function to solve<a name="line.239"></a>
243 <FONT color="green">240</FONT> * @param x0 old approximation for the root<a name="line.240"></a>
244 <FONT color="green">241</FONT> * @param y0 function value at the approximation for the root<a name="line.241"></a>
245 <FONT color="green">242</FONT> * @param x1 last calculated approximation for the root<a name="line.242"></a>
246 <FONT color="green">243</FONT> * @param y1 function value at the last calculated approximation<a name="line.243"></a>
247 <FONT color="green">244</FONT> * for the root<a name="line.244"></a>
248 <FONT color="green">245</FONT> * @param x2 bracket point (must be set to x0 if no bracket point is<a name="line.245"></a>
249 <FONT color="green">246</FONT> * known, this will force starting with linear interpolation)<a name="line.246"></a>
250 <FONT color="green">247</FONT> * @param y2 function value at the bracket point.<a name="line.247"></a>
251 <FONT color="green">248</FONT> * @return the value where the function is zero<a name="line.248"></a>
252 <FONT color="green">249</FONT> * @throws MaxIterationsExceededException if the maximum iteration count<a name="line.249"></a>
253 <FONT color="green">250</FONT> * is exceeded<a name="line.250"></a>
254 <FONT color="green">251</FONT> * @throws FunctionEvaluationException if an error occurs evaluating<a name="line.251"></a>
255 <FONT color="green">252</FONT> * the function<a name="line.252"></a>
256 <FONT color="green">253</FONT> */<a name="line.253"></a>
257 <FONT color="green">254</FONT> private double solve(final UnivariateRealFunction f,<a name="line.254"></a>
258 <FONT color="green">255</FONT> double x0, double y0,<a name="line.255"></a>
259 <FONT color="green">256</FONT> double x1, double y1,<a name="line.256"></a>
260 <FONT color="green">257</FONT> double x2, double y2)<a name="line.257"></a>
261 <FONT color="green">258</FONT> throws MaxIterationsExceededException, FunctionEvaluationException {<a name="line.258"></a>
262 <FONT color="green">259</FONT> <a name="line.259"></a>
263 <FONT color="green">260</FONT> double delta = x1 - x0;<a name="line.260"></a>
264 <FONT color="green">261</FONT> double oldDelta = delta;<a name="line.261"></a>
265 <FONT color="green">262</FONT> <a name="line.262"></a>
266 <FONT color="green">263</FONT> int i = 0;<a name="line.263"></a>
267 <FONT color="green">264</FONT> while (i &lt; maximalIterationCount) {<a name="line.264"></a>
268 <FONT color="green">265</FONT> if (Math.abs(y2) &lt; Math.abs(y1)) {<a name="line.265"></a>
269 <FONT color="green">266</FONT> // use the bracket point if is better than last approximation<a name="line.266"></a>
270 <FONT color="green">267</FONT> x0 = x1;<a name="line.267"></a>
271 <FONT color="green">268</FONT> x1 = x2;<a name="line.268"></a>
272 <FONT color="green">269</FONT> x2 = x0;<a name="line.269"></a>
273 <FONT color="green">270</FONT> y0 = y1;<a name="line.270"></a>
274 <FONT color="green">271</FONT> y1 = y2;<a name="line.271"></a>
275 <FONT color="green">272</FONT> y2 = y0;<a name="line.272"></a>
276 <FONT color="green">273</FONT> }<a name="line.273"></a>
277 <FONT color="green">274</FONT> if (Math.abs(y1) &lt;= functionValueAccuracy) {<a name="line.274"></a>
278 <FONT color="green">275</FONT> // Avoid division by very small values. Assume<a name="line.275"></a>
279 <FONT color="green">276</FONT> // the iteration has converged (the problem may<a name="line.276"></a>
280 <FONT color="green">277</FONT> // still be ill conditioned)<a name="line.277"></a>
281 <FONT color="green">278</FONT> setResult(x1, i);<a name="line.278"></a>
282 <FONT color="green">279</FONT> return result;<a name="line.279"></a>
283 <FONT color="green">280</FONT> }<a name="line.280"></a>
284 <FONT color="green">281</FONT> double dx = x2 - x1;<a name="line.281"></a>
285 <FONT color="green">282</FONT> double tolerance =<a name="line.282"></a>
286 <FONT color="green">283</FONT> Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy);<a name="line.283"></a>
287 <FONT color="green">284</FONT> if (Math.abs(dx) &lt;= tolerance) {<a name="line.284"></a>
288 <FONT color="green">285</FONT> setResult(x1, i);<a name="line.285"></a>
289 <FONT color="green">286</FONT> return result;<a name="line.286"></a>
290 <FONT color="green">287</FONT> }<a name="line.287"></a>
291 <FONT color="green">288</FONT> if ((Math.abs(oldDelta) &lt; tolerance) ||<a name="line.288"></a>
292 <FONT color="green">289</FONT> (Math.abs(y0) &lt;= Math.abs(y1))) {<a name="line.289"></a>
293 <FONT color="green">290</FONT> // Force bisection.<a name="line.290"></a>
294 <FONT color="green">291</FONT> delta = 0.5 * dx;<a name="line.291"></a>
295 <FONT color="green">292</FONT> oldDelta = delta;<a name="line.292"></a>
296 <FONT color="green">293</FONT> } else {<a name="line.293"></a>
297 <FONT color="green">294</FONT> double r3 = y1 / y0;<a name="line.294"></a>
298 <FONT color="green">295</FONT> double p;<a name="line.295"></a>
299 <FONT color="green">296</FONT> double p1;<a name="line.296"></a>
300 <FONT color="green">297</FONT> // the equality test (x0 == x2) is intentional,<a name="line.297"></a>
301 <FONT color="green">298</FONT> // it is part of the original Brent's method,<a name="line.298"></a>
302 <FONT color="green">299</FONT> // it should NOT be replaced by proximity test<a name="line.299"></a>
303 <FONT color="green">300</FONT> if (x0 == x2) {<a name="line.300"></a>
304 <FONT color="green">301</FONT> // Linear interpolation.<a name="line.301"></a>
305 <FONT color="green">302</FONT> p = dx * r3;<a name="line.302"></a>
306 <FONT color="green">303</FONT> p1 = 1.0 - r3;<a name="line.303"></a>
307 <FONT color="green">304</FONT> } else {<a name="line.304"></a>
308 <FONT color="green">305</FONT> // Inverse quadratic interpolation.<a name="line.305"></a>
309 <FONT color="green">306</FONT> double r1 = y0 / y2;<a name="line.306"></a>
310 <FONT color="green">307</FONT> double r2 = y1 / y2;<a name="line.307"></a>
311 <FONT color="green">308</FONT> p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));<a name="line.308"></a>
312 <FONT color="green">309</FONT> p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);<a name="line.309"></a>
313 <FONT color="green">310</FONT> }<a name="line.310"></a>
314 <FONT color="green">311</FONT> if (p &gt; 0.0) {<a name="line.311"></a>
315 <FONT color="green">312</FONT> p1 = -p1;<a name="line.312"></a>
316 <FONT color="green">313</FONT> } else {<a name="line.313"></a>
317 <FONT color="green">314</FONT> p = -p;<a name="line.314"></a>
318 <FONT color="green">315</FONT> }<a name="line.315"></a>
319 <FONT color="green">316</FONT> if (2.0 * p &gt;= 1.5 * dx * p1 - Math.abs(tolerance * p1) ||<a name="line.316"></a>
320 <FONT color="green">317</FONT> p &gt;= Math.abs(0.5 * oldDelta * p1)) {<a name="line.317"></a>
321 <FONT color="green">318</FONT> // Inverse quadratic interpolation gives a value<a name="line.318"></a>
322 <FONT color="green">319</FONT> // in the wrong direction, or progress is slow.<a name="line.319"></a>
323 <FONT color="green">320</FONT> // Fall back to bisection.<a name="line.320"></a>
324 <FONT color="green">321</FONT> delta = 0.5 * dx;<a name="line.321"></a>
325 <FONT color="green">322</FONT> oldDelta = delta;<a name="line.322"></a>
326 <FONT color="green">323</FONT> } else {<a name="line.323"></a>
327 <FONT color="green">324</FONT> oldDelta = delta;<a name="line.324"></a>
328 <FONT color="green">325</FONT> delta = p / p1;<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> // Save old X1, Y1<a name="line.328"></a>
332 <FONT color="green">329</FONT> x0 = x1;<a name="line.329"></a>
333 <FONT color="green">330</FONT> y0 = y1;<a name="line.330"></a>
334 <FONT color="green">331</FONT> // Compute new X1, Y1<a name="line.331"></a>
335 <FONT color="green">332</FONT> if (Math.abs(delta) &gt; tolerance) {<a name="line.332"></a>
336 <FONT color="green">333</FONT> x1 = x1 + delta;<a name="line.333"></a>
337 <FONT color="green">334</FONT> } else if (dx &gt; 0.0) {<a name="line.334"></a>
338 <FONT color="green">335</FONT> x1 = x1 + 0.5 * tolerance;<a name="line.335"></a>
339 <FONT color="green">336</FONT> } else if (dx &lt;= 0.0) {<a name="line.336"></a>
340 <FONT color="green">337</FONT> x1 = x1 - 0.5 * tolerance;<a name="line.337"></a>
341 <FONT color="green">338</FONT> }<a name="line.338"></a>
342 <FONT color="green">339</FONT> y1 = f.value(x1);<a name="line.339"></a>
343 <FONT color="green">340</FONT> if ((y1 &gt; 0) == (y2 &gt; 0)) {<a name="line.340"></a>
344 <FONT color="green">341</FONT> x2 = x0;<a name="line.341"></a>
345 <FONT color="green">342</FONT> y2 = y0;<a name="line.342"></a>
346 <FONT color="green">343</FONT> delta = x1 - x0;<a name="line.343"></a>
347 <FONT color="green">344</FONT> oldDelta = delta;<a name="line.344"></a>
348 <FONT color="green">345</FONT> }<a name="line.345"></a>
349 <FONT color="green">346</FONT> i++;<a name="line.346"></a>
350 <FONT color="green">347</FONT> }<a name="line.347"></a>
351 <FONT color="green">348</FONT> throw new MaxIterationsExceededException(maximalIterationCount);<a name="line.348"></a>
352 <FONT color="green">349</FONT> }<a name="line.349"></a>
353 <FONT color="green">350</FONT> }<a name="line.350"></a>
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414 </PRE>
415 </BODY>
416 </HTML>