comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/analysis/solvers/RiddersSolver.html @ 13:cbf34dd4d7e6

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> import org.apache.commons.math.ConvergenceException;<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.MaxIterationsExceededException;<a name="line.21"></a>
25 <FONT color="green">022</FONT> import org.apache.commons.math.analysis.UnivariateRealFunction;<a name="line.22"></a>
26 <FONT color="green">023</FONT> import org.apache.commons.math.util.MathUtils;<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/RiddersMethod.html"&gt;<a name="line.26"></a>
30 <FONT color="green">027</FONT> * Ridders' Method&lt;/a&gt; for root finding of real univariate functions. For<a name="line.27"></a>
31 <FONT color="green">028</FONT> * reference, see C. Ridders, &lt;i&gt;A new algorithm for computing a single root<a name="line.28"></a>
32 <FONT color="green">029</FONT> * of a real continuous function &lt;/i&gt;, IEEE Transactions on Circuits and<a name="line.29"></a>
33 <FONT color="green">030</FONT> * Systems, 26 (1979), 979 - 980.<a name="line.30"></a>
34 <FONT color="green">031</FONT> * &lt;p&gt;<a name="line.31"></a>
35 <FONT color="green">032</FONT> * The function should be continuous but not necessarily smooth.&lt;/p&gt;<a name="line.32"></a>
36 <FONT color="green">033</FONT> *<a name="line.33"></a>
37 <FONT color="green">034</FONT> * @version $Revision: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 2009) $<a name="line.34"></a>
38 <FONT color="green">035</FONT> * @since 1.2<a name="line.35"></a>
39 <FONT color="green">036</FONT> */<a name="line.36"></a>
40 <FONT color="green">037</FONT> public class RiddersSolver extends UnivariateRealSolverImpl {<a name="line.37"></a>
41 <FONT color="green">038</FONT> <a name="line.38"></a>
42 <FONT color="green">039</FONT> /**<a name="line.39"></a>
43 <FONT color="green">040</FONT> * Construct a solver for the given function.<a name="line.40"></a>
44 <FONT color="green">041</FONT> *<a name="line.41"></a>
45 <FONT color="green">042</FONT> * @param f function to solve<a name="line.42"></a>
46 <FONT color="green">043</FONT> * @deprecated as of 2.0 the function to solve is passed as an argument<a name="line.43"></a>
47 <FONT color="green">044</FONT> * to the {@link #solve(UnivariateRealFunction, double, double)} or<a name="line.44"></a>
48 <FONT color="green">045</FONT> * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}<a name="line.45"></a>
49 <FONT color="green">046</FONT> * method.<a name="line.46"></a>
50 <FONT color="green">047</FONT> */<a name="line.47"></a>
51 <FONT color="green">048</FONT> @Deprecated<a name="line.48"></a>
52 <FONT color="green">049</FONT> public RiddersSolver(UnivariateRealFunction f) {<a name="line.49"></a>
53 <FONT color="green">050</FONT> super(f, 100, 1E-6);<a name="line.50"></a>
54 <FONT color="green">051</FONT> }<a name="line.51"></a>
55 <FONT color="green">052</FONT> <a name="line.52"></a>
56 <FONT color="green">053</FONT> /**<a name="line.53"></a>
57 <FONT color="green">054</FONT> * Construct a solver.<a name="line.54"></a>
58 <FONT color="green">055</FONT> */<a name="line.55"></a>
59 <FONT color="green">056</FONT> public RiddersSolver() {<a name="line.56"></a>
60 <FONT color="green">057</FONT> super(100, 1E-6);<a name="line.57"></a>
61 <FONT color="green">058</FONT> }<a name="line.58"></a>
62 <FONT color="green">059</FONT> <a name="line.59"></a>
63 <FONT color="green">060</FONT> /** {@inheritDoc} */<a name="line.60"></a>
64 <FONT color="green">061</FONT> @Deprecated<a name="line.61"></a>
65 <FONT color="green">062</FONT> public double solve(final double min, final double max)<a name="line.62"></a>
66 <FONT color="green">063</FONT> throws ConvergenceException, FunctionEvaluationException {<a name="line.63"></a>
67 <FONT color="green">064</FONT> return solve(f, min, max);<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> /** {@inheritDoc} */<a name="line.67"></a>
71 <FONT color="green">068</FONT> @Deprecated<a name="line.68"></a>
72 <FONT color="green">069</FONT> public double solve(final double min, final double max, final double initial)<a name="line.69"></a>
73 <FONT color="green">070</FONT> throws ConvergenceException, FunctionEvaluationException {<a name="line.70"></a>
74 <FONT color="green">071</FONT> return solve(f, min, max, initial);<a name="line.71"></a>
75 <FONT color="green">072</FONT> }<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> * Find a root in the given interval with initial value.<a name="line.75"></a>
79 <FONT color="green">076</FONT> * &lt;p&gt;<a name="line.76"></a>
80 <FONT color="green">077</FONT> * Requires bracketing condition.&lt;/p&gt;<a name="line.77"></a>
81 <FONT color="green">078</FONT> *<a name="line.78"></a>
82 <FONT color="green">079</FONT> * @param f the function to solve<a name="line.79"></a>
83 <FONT color="green">080</FONT> * @param min the lower bound for the interval<a name="line.80"></a>
84 <FONT color="green">081</FONT> * @param max the upper bound for the interval<a name="line.81"></a>
85 <FONT color="green">082</FONT> * @param initial the start value to use<a name="line.82"></a>
86 <FONT color="green">083</FONT> * @return the point at which the function value is zero<a name="line.83"></a>
87 <FONT color="green">084</FONT> * @throws MaxIterationsExceededException if the maximum iteration count is exceeded<a name="line.84"></a>
88 <FONT color="green">085</FONT> * @throws FunctionEvaluationException if an error occurs evaluating the<a name="line.85"></a>
89 <FONT color="green">086</FONT> * function<a name="line.86"></a>
90 <FONT color="green">087</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.87"></a>
91 <FONT color="green">088</FONT> */<a name="line.88"></a>
92 <FONT color="green">089</FONT> public double solve(final UnivariateRealFunction f,<a name="line.89"></a>
93 <FONT color="green">090</FONT> final double min, final double max, final double initial)<a name="line.90"></a>
94 <FONT color="green">091</FONT> throws MaxIterationsExceededException, FunctionEvaluationException {<a name="line.91"></a>
95 <FONT color="green">092</FONT> <a name="line.92"></a>
96 <FONT color="green">093</FONT> // check for zeros before verifying bracketing<a name="line.93"></a>
97 <FONT color="green">094</FONT> if (f.value(min) == 0.0) { return min; }<a name="line.94"></a>
98 <FONT color="green">095</FONT> if (f.value(max) == 0.0) { return max; }<a name="line.95"></a>
99 <FONT color="green">096</FONT> if (f.value(initial) == 0.0) { return initial; }<a name="line.96"></a>
100 <FONT color="green">097</FONT> <a name="line.97"></a>
101 <FONT color="green">098</FONT> verifyBracketing(min, max, f);<a name="line.98"></a>
102 <FONT color="green">099</FONT> verifySequence(min, initial, max);<a name="line.99"></a>
103 <FONT color="green">100</FONT> if (isBracketing(min, initial, f)) {<a name="line.100"></a>
104 <FONT color="green">101</FONT> return solve(f, min, initial);<a name="line.101"></a>
105 <FONT color="green">102</FONT> } else {<a name="line.102"></a>
106 <FONT color="green">103</FONT> return solve(f, initial, max);<a name="line.103"></a>
107 <FONT color="green">104</FONT> }<a name="line.104"></a>
108 <FONT color="green">105</FONT> }<a name="line.105"></a>
109 <FONT color="green">106</FONT> <a name="line.106"></a>
110 <FONT color="green">107</FONT> /**<a name="line.107"></a>
111 <FONT color="green">108</FONT> * Find a root in the given interval.<a name="line.108"></a>
112 <FONT color="green">109</FONT> * &lt;p&gt;<a name="line.109"></a>
113 <FONT color="green">110</FONT> * Requires bracketing condition.&lt;/p&gt;<a name="line.110"></a>
114 <FONT color="green">111</FONT> *<a name="line.111"></a>
115 <FONT color="green">112</FONT> * @param f the function to solve<a name="line.112"></a>
116 <FONT color="green">113</FONT> * @param min the lower bound for the interval<a name="line.113"></a>
117 <FONT color="green">114</FONT> * @param max the upper bound for the interval<a name="line.114"></a>
118 <FONT color="green">115</FONT> * @return the point at which the function value is zero<a name="line.115"></a>
119 <FONT color="green">116</FONT> * @throws MaxIterationsExceededException if the maximum iteration count is exceeded<a name="line.116"></a>
120 <FONT color="green">117</FONT> * @throws FunctionEvaluationException if an error occurs evaluating the<a name="line.117"></a>
121 <FONT color="green">118</FONT> * function<a name="line.118"></a>
122 <FONT color="green">119</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.119"></a>
123 <FONT color="green">120</FONT> */<a name="line.120"></a>
124 <FONT color="green">121</FONT> public double solve(final UnivariateRealFunction f,<a name="line.121"></a>
125 <FONT color="green">122</FONT> final double min, final double max)<a name="line.122"></a>
126 <FONT color="green">123</FONT> throws MaxIterationsExceededException, FunctionEvaluationException {<a name="line.123"></a>
127 <FONT color="green">124</FONT> <a name="line.124"></a>
128 <FONT color="green">125</FONT> // [x1, x2] is the bracketing interval in each iteration<a name="line.125"></a>
129 <FONT color="green">126</FONT> // x3 is the midpoint of [x1, x2]<a name="line.126"></a>
130 <FONT color="green">127</FONT> // x is the new root approximation and an endpoint of the new interval<a name="line.127"></a>
131 <FONT color="green">128</FONT> double x1 = min;<a name="line.128"></a>
132 <FONT color="green">129</FONT> double y1 = f.value(x1);<a name="line.129"></a>
133 <FONT color="green">130</FONT> double x2 = max;<a name="line.130"></a>
134 <FONT color="green">131</FONT> double y2 = f.value(x2);<a name="line.131"></a>
135 <FONT color="green">132</FONT> <a name="line.132"></a>
136 <FONT color="green">133</FONT> // check for zeros before verifying bracketing<a name="line.133"></a>
137 <FONT color="green">134</FONT> if (y1 == 0.0) {<a name="line.134"></a>
138 <FONT color="green">135</FONT> return min;<a name="line.135"></a>
139 <FONT color="green">136</FONT> }<a name="line.136"></a>
140 <FONT color="green">137</FONT> if (y2 == 0.0) {<a name="line.137"></a>
141 <FONT color="green">138</FONT> return max;<a name="line.138"></a>
142 <FONT color="green">139</FONT> }<a name="line.139"></a>
143 <FONT color="green">140</FONT> verifyBracketing(min, max, f);<a name="line.140"></a>
144 <FONT color="green">141</FONT> <a name="line.141"></a>
145 <FONT color="green">142</FONT> int i = 1;<a name="line.142"></a>
146 <FONT color="green">143</FONT> double oldx = Double.POSITIVE_INFINITY;<a name="line.143"></a>
147 <FONT color="green">144</FONT> while (i &lt;= maximalIterationCount) {<a name="line.144"></a>
148 <FONT color="green">145</FONT> // calculate the new root approximation<a name="line.145"></a>
149 <FONT color="green">146</FONT> final double x3 = 0.5 * (x1 + x2);<a name="line.146"></a>
150 <FONT color="green">147</FONT> final double y3 = f.value(x3);<a name="line.147"></a>
151 <FONT color="green">148</FONT> if (Math.abs(y3) &lt;= functionValueAccuracy) {<a name="line.148"></a>
152 <FONT color="green">149</FONT> setResult(x3, i);<a name="line.149"></a>
153 <FONT color="green">150</FONT> return result;<a name="line.150"></a>
154 <FONT color="green">151</FONT> }<a name="line.151"></a>
155 <FONT color="green">152</FONT> final double delta = 1 - (y1 * y2) / (y3 * y3); // delta &gt; 1 due to bracketing<a name="line.152"></a>
156 <FONT color="green">153</FONT> final double correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *<a name="line.153"></a>
157 <FONT color="green">154</FONT> (x3 - x1) / Math.sqrt(delta);<a name="line.154"></a>
158 <FONT color="green">155</FONT> final double x = x3 - correction; // correction != 0<a name="line.155"></a>
159 <FONT color="green">156</FONT> final double y = f.value(x);<a name="line.156"></a>
160 <FONT color="green">157</FONT> <a name="line.157"></a>
161 <FONT color="green">158</FONT> // check for convergence<a name="line.158"></a>
162 <FONT color="green">159</FONT> final double tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);<a name="line.159"></a>
163 <FONT color="green">160</FONT> if (Math.abs(x - oldx) &lt;= tolerance) {<a name="line.160"></a>
164 <FONT color="green">161</FONT> setResult(x, i);<a name="line.161"></a>
165 <FONT color="green">162</FONT> return result;<a name="line.162"></a>
166 <FONT color="green">163</FONT> }<a name="line.163"></a>
167 <FONT color="green">164</FONT> if (Math.abs(y) &lt;= functionValueAccuracy) {<a name="line.164"></a>
168 <FONT color="green">165</FONT> setResult(x, i);<a name="line.165"></a>
169 <FONT color="green">166</FONT> return result;<a name="line.166"></a>
170 <FONT color="green">167</FONT> }<a name="line.167"></a>
171 <FONT color="green">168</FONT> <a name="line.168"></a>
172 <FONT color="green">169</FONT> // prepare the new interval for next iteration<a name="line.169"></a>
173 <FONT color="green">170</FONT> // Ridders' method guarantees x1 &lt; x &lt; x2<a name="line.170"></a>
174 <FONT color="green">171</FONT> if (correction &gt; 0.0) { // x1 &lt; x &lt; x3<a name="line.171"></a>
175 <FONT color="green">172</FONT> if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) {<a name="line.172"></a>
176 <FONT color="green">173</FONT> x2 = x;<a name="line.173"></a>
177 <FONT color="green">174</FONT> y2 = y;<a name="line.174"></a>
178 <FONT color="green">175</FONT> } else {<a name="line.175"></a>
179 <FONT color="green">176</FONT> x1 = x;<a name="line.176"></a>
180 <FONT color="green">177</FONT> x2 = x3;<a name="line.177"></a>
181 <FONT color="green">178</FONT> y1 = y;<a name="line.178"></a>
182 <FONT color="green">179</FONT> y2 = y3;<a name="line.179"></a>
183 <FONT color="green">180</FONT> }<a name="line.180"></a>
184 <FONT color="green">181</FONT> } else { // x3 &lt; x &lt; x2<a name="line.181"></a>
185 <FONT color="green">182</FONT> if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) {<a name="line.182"></a>
186 <FONT color="green">183</FONT> x1 = x;<a name="line.183"></a>
187 <FONT color="green">184</FONT> y1 = y;<a name="line.184"></a>
188 <FONT color="green">185</FONT> } else {<a name="line.185"></a>
189 <FONT color="green">186</FONT> x1 = x3;<a name="line.186"></a>
190 <FONT color="green">187</FONT> x2 = x;<a name="line.187"></a>
191 <FONT color="green">188</FONT> y1 = y3;<a name="line.188"></a>
192 <FONT color="green">189</FONT> y2 = y;<a name="line.189"></a>
193 <FONT color="green">190</FONT> }<a name="line.190"></a>
194 <FONT color="green">191</FONT> }<a name="line.191"></a>
195 <FONT color="green">192</FONT> oldx = x;<a name="line.192"></a>
196 <FONT color="green">193</FONT> i++;<a name="line.193"></a>
197 <FONT color="green">194</FONT> }<a name="line.194"></a>
198 <FONT color="green">195</FONT> throw new MaxIterationsExceededException(maximalIterationCount);<a name="line.195"></a>
199 <FONT color="green">196</FONT> }<a name="line.196"></a>
200 <FONT color="green">197</FONT> }<a name="line.197"></a>
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261 </PRE>
262 </BODY>
263 </HTML>