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

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/analysis/solvers/LaguerreSolver.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> 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.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> import org.apache.commons.math.analysis.polynomials.PolynomialFunction;<a name="line.24"></a>
28 <FONT color="green">025</FONT> import org.apache.commons.math.complex.Complex;<a name="line.25"></a>
29 <FONT color="green">026</FONT> <a name="line.26"></a>
30 <FONT color="green">027</FONT> /**<a name="line.27"></a>
31 <FONT color="green">028</FONT> * Implements the &lt;a href="http://mathworld.wolfram.com/LaguerresMethod.html"&gt;<a name="line.28"></a>
32 <FONT color="green">029</FONT> * Laguerre's Method&lt;/a&gt; for root finding of real coefficient polynomials.<a name="line.29"></a>
33 <FONT color="green">030</FONT> * For reference, see &lt;b&gt;A First Course in Numerical Analysis&lt;/b&gt;,<a name="line.30"></a>
34 <FONT color="green">031</FONT> * ISBN 048641454X, chapter 8.<a name="line.31"></a>
35 <FONT color="green">032</FONT> * &lt;p&gt;<a name="line.32"></a>
36 <FONT color="green">033</FONT> * Laguerre's method is global in the sense that it can start with any initial<a name="line.33"></a>
37 <FONT color="green">034</FONT> * approximation and be able to solve all roots from that point.&lt;/p&gt;<a name="line.34"></a>
38 <FONT color="green">035</FONT> *<a name="line.35"></a>
39 <FONT color="green">036</FONT> * @version $Revision: 922708 $ $Date: 2010-03-13 20:15:47 -0500 (Sat, 13 Mar 2010) $<a name="line.36"></a>
40 <FONT color="green">037</FONT> * @since 1.2<a name="line.37"></a>
41 <FONT color="green">038</FONT> */<a name="line.38"></a>
42 <FONT color="green">039</FONT> public class LaguerreSolver extends UnivariateRealSolverImpl {<a name="line.39"></a>
43 <FONT color="green">040</FONT> <a name="line.40"></a>
44 <FONT color="green">041</FONT> /** Message for non-polynomial function. */<a name="line.41"></a>
45 <FONT color="green">042</FONT> private static final String NON_POLYNOMIAL_FUNCTION_MESSAGE =<a name="line.42"></a>
46 <FONT color="green">043</FONT> "function is not polynomial";<a name="line.43"></a>
47 <FONT color="green">044</FONT> <a name="line.44"></a>
48 <FONT color="green">045</FONT> /** Message for non-positive degree. */<a name="line.45"></a>
49 <FONT color="green">046</FONT> private static final String NON_POSITIVE_DEGREE_MESSAGE =<a name="line.46"></a>
50 <FONT color="green">047</FONT> "polynomial degree must be positive: degree={0}";<a name="line.47"></a>
51 <FONT color="green">048</FONT> <a name="line.48"></a>
52 <FONT color="green">049</FONT> /** polynomial function to solve.<a name="line.49"></a>
53 <FONT color="green">050</FONT> * @deprecated as of 2.0 the function is not stored anymore in the instance<a name="line.50"></a>
54 <FONT color="green">051</FONT> */<a name="line.51"></a>
55 <FONT color="green">052</FONT> @Deprecated<a name="line.52"></a>
56 <FONT color="green">053</FONT> private final PolynomialFunction p;<a name="line.53"></a>
57 <FONT color="green">054</FONT> <a name="line.54"></a>
58 <FONT color="green">055</FONT> /**<a name="line.55"></a>
59 <FONT color="green">056</FONT> * Construct a solver for the given function.<a name="line.56"></a>
60 <FONT color="green">057</FONT> *<a name="line.57"></a>
61 <FONT color="green">058</FONT> * @param f function to solve<a name="line.58"></a>
62 <FONT color="green">059</FONT> * @throws IllegalArgumentException if function is not polynomial<a name="line.59"></a>
63 <FONT color="green">060</FONT> * @deprecated as of 2.0 the function to solve is passed as an argument<a name="line.60"></a>
64 <FONT color="green">061</FONT> * to the {@link #solve(UnivariateRealFunction, double, double)} or<a name="line.61"></a>
65 <FONT color="green">062</FONT> * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}<a name="line.62"></a>
66 <FONT color="green">063</FONT> * method.<a name="line.63"></a>
67 <FONT color="green">064</FONT> */<a name="line.64"></a>
68 <FONT color="green">065</FONT> @Deprecated<a name="line.65"></a>
69 <FONT color="green">066</FONT> public LaguerreSolver(UnivariateRealFunction f) throws<a name="line.66"></a>
70 <FONT color="green">067</FONT> IllegalArgumentException {<a name="line.67"></a>
71 <FONT color="green">068</FONT> super(f, 100, 1E-6);<a name="line.68"></a>
72 <FONT color="green">069</FONT> if (f instanceof PolynomialFunction) {<a name="line.69"></a>
73 <FONT color="green">070</FONT> p = (PolynomialFunction) f;<a name="line.70"></a>
74 <FONT color="green">071</FONT> } else {<a name="line.71"></a>
75 <FONT color="green">072</FONT> throw MathRuntimeException.createIllegalArgumentException(NON_POLYNOMIAL_FUNCTION_MESSAGE);<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> /**<a name="line.76"></a>
80 <FONT color="green">077</FONT> * Construct a solver.<a name="line.77"></a>
81 <FONT color="green">078</FONT> */<a name="line.78"></a>
82 <FONT color="green">079</FONT> public LaguerreSolver() {<a name="line.79"></a>
83 <FONT color="green">080</FONT> super(100, 1E-6);<a name="line.80"></a>
84 <FONT color="green">081</FONT> p = null;<a name="line.81"></a>
85 <FONT color="green">082</FONT> }<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> * Returns a copy of the polynomial function.<a name="line.85"></a>
89 <FONT color="green">086</FONT> *<a name="line.86"></a>
90 <FONT color="green">087</FONT> * @return a fresh copy of the polynomial function<a name="line.87"></a>
91 <FONT color="green">088</FONT> * @deprecated as of 2.0 the function is not stored anymore within the instance.<a name="line.88"></a>
92 <FONT color="green">089</FONT> */<a name="line.89"></a>
93 <FONT color="green">090</FONT> @Deprecated<a name="line.90"></a>
94 <FONT color="green">091</FONT> public PolynomialFunction getPolynomialFunction() {<a name="line.91"></a>
95 <FONT color="green">092</FONT> return new PolynomialFunction(p.getCoefficients());<a name="line.92"></a>
96 <FONT color="green">093</FONT> }<a name="line.93"></a>
97 <FONT color="green">094</FONT> <a name="line.94"></a>
98 <FONT color="green">095</FONT> /** {@inheritDoc} */<a name="line.95"></a>
99 <FONT color="green">096</FONT> @Deprecated<a name="line.96"></a>
100 <FONT color="green">097</FONT> public double solve(final double min, final double max)<a name="line.97"></a>
101 <FONT color="green">098</FONT> throws ConvergenceException, FunctionEvaluationException {<a name="line.98"></a>
102 <FONT color="green">099</FONT> return solve(p, min, max);<a name="line.99"></a>
103 <FONT color="green">100</FONT> }<a name="line.100"></a>
104 <FONT color="green">101</FONT> <a name="line.101"></a>
105 <FONT color="green">102</FONT> /** {@inheritDoc} */<a name="line.102"></a>
106 <FONT color="green">103</FONT> @Deprecated<a name="line.103"></a>
107 <FONT color="green">104</FONT> public double solve(final double min, final double max, final double initial)<a name="line.104"></a>
108 <FONT color="green">105</FONT> throws ConvergenceException, FunctionEvaluationException {<a name="line.105"></a>
109 <FONT color="green">106</FONT> return solve(p, min, max, initial);<a name="line.106"></a>
110 <FONT color="green">107</FONT> }<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> * Find a real root in the given interval with initial value.<a name="line.110"></a>
114 <FONT color="green">111</FONT> * &lt;p&gt;<a name="line.111"></a>
115 <FONT color="green">112</FONT> * Requires bracketing condition.&lt;/p&gt;<a name="line.112"></a>
116 <FONT color="green">113</FONT> *<a name="line.113"></a>
117 <FONT color="green">114</FONT> * @param f function to solve (must be polynomial)<a name="line.114"></a>
118 <FONT color="green">115</FONT> * @param min the lower bound for the interval<a name="line.115"></a>
119 <FONT color="green">116</FONT> * @param max the upper bound for the interval<a name="line.116"></a>
120 <FONT color="green">117</FONT> * @param initial the start value to use<a name="line.117"></a>
121 <FONT color="green">118</FONT> * @return the point at which the function value is zero<a name="line.118"></a>
122 <FONT color="green">119</FONT> * @throws ConvergenceException if the maximum iteration count is exceeded<a name="line.119"></a>
123 <FONT color="green">120</FONT> * or the solver detects convergence problems otherwise<a name="line.120"></a>
124 <FONT color="green">121</FONT> * @throws FunctionEvaluationException if an error occurs evaluating the<a name="line.121"></a>
125 <FONT color="green">122</FONT> * function<a name="line.122"></a>
126 <FONT color="green">123</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.123"></a>
127 <FONT color="green">124</FONT> */<a name="line.124"></a>
128 <FONT color="green">125</FONT> public double solve(final UnivariateRealFunction f,<a name="line.125"></a>
129 <FONT color="green">126</FONT> final double min, final double max, final double initial)<a name="line.126"></a>
130 <FONT color="green">127</FONT> throws ConvergenceException, FunctionEvaluationException {<a name="line.127"></a>
131 <FONT color="green">128</FONT> <a name="line.128"></a>
132 <FONT color="green">129</FONT> // check for zeros before verifying bracketing<a name="line.129"></a>
133 <FONT color="green">130</FONT> if (f.value(min) == 0.0) {<a name="line.130"></a>
134 <FONT color="green">131</FONT> return min;<a name="line.131"></a>
135 <FONT color="green">132</FONT> }<a name="line.132"></a>
136 <FONT color="green">133</FONT> if (f.value(max) == 0.0) {<a name="line.133"></a>
137 <FONT color="green">134</FONT> return max;<a name="line.134"></a>
138 <FONT color="green">135</FONT> }<a name="line.135"></a>
139 <FONT color="green">136</FONT> if (f.value(initial) == 0.0) {<a name="line.136"></a>
140 <FONT color="green">137</FONT> return initial;<a name="line.137"></a>
141 <FONT color="green">138</FONT> }<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> verifySequence(min, initial, max);<a name="line.141"></a>
145 <FONT color="green">142</FONT> if (isBracketing(min, initial, f)) {<a name="line.142"></a>
146 <FONT color="green">143</FONT> return solve(f, min, initial);<a name="line.143"></a>
147 <FONT color="green">144</FONT> } else {<a name="line.144"></a>
148 <FONT color="green">145</FONT> return solve(f, initial, max);<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> }<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> * Find a real root in the given interval.<a name="line.151"></a>
155 <FONT color="green">152</FONT> * &lt;p&gt;<a name="line.152"></a>
156 <FONT color="green">153</FONT> * Despite the bracketing condition, the root returned by solve(Complex[],<a name="line.153"></a>
157 <FONT color="green">154</FONT> * Complex) may not be a real zero inside [min, max]. For example,<a name="line.154"></a>
158 <FONT color="green">155</FONT> * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try<a name="line.155"></a>
159 <FONT color="green">156</FONT> * another initial value, or, as we did here, call solveAll() to obtain<a name="line.156"></a>
160 <FONT color="green">157</FONT> * all roots and pick up the one that we're looking for.&lt;/p&gt;<a name="line.157"></a>
161 <FONT color="green">158</FONT> *<a name="line.158"></a>
162 <FONT color="green">159</FONT> * @param f the function to solve<a name="line.159"></a>
163 <FONT color="green">160</FONT> * @param min the lower bound for the interval<a name="line.160"></a>
164 <FONT color="green">161</FONT> * @param max the upper bound for the interval<a name="line.161"></a>
165 <FONT color="green">162</FONT> * @return the point at which the function value is zero<a name="line.162"></a>
166 <FONT color="green">163</FONT> * @throws ConvergenceException if the maximum iteration count is exceeded<a name="line.163"></a>
167 <FONT color="green">164</FONT> * or the solver detects convergence problems otherwise<a name="line.164"></a>
168 <FONT color="green">165</FONT> * @throws FunctionEvaluationException if an error occurs evaluating the<a name="line.165"></a>
169 <FONT color="green">166</FONT> * function<a name="line.166"></a>
170 <FONT color="green">167</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.167"></a>
171 <FONT color="green">168</FONT> */<a name="line.168"></a>
172 <FONT color="green">169</FONT> public double solve(final UnivariateRealFunction f,<a name="line.169"></a>
173 <FONT color="green">170</FONT> final double min, final double max)<a name="line.170"></a>
174 <FONT color="green">171</FONT> throws ConvergenceException, FunctionEvaluationException {<a name="line.171"></a>
175 <FONT color="green">172</FONT> <a name="line.172"></a>
176 <FONT color="green">173</FONT> // check function type<a name="line.173"></a>
177 <FONT color="green">174</FONT> if (!(f instanceof PolynomialFunction)) {<a name="line.174"></a>
178 <FONT color="green">175</FONT> throw MathRuntimeException.createIllegalArgumentException(NON_POLYNOMIAL_FUNCTION_MESSAGE);<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> // check for zeros before verifying bracketing<a name="line.178"></a>
182 <FONT color="green">179</FONT> if (f.value(min) == 0.0) { return min; }<a name="line.179"></a>
183 <FONT color="green">180</FONT> if (f.value(max) == 0.0) { return max; }<a name="line.180"></a>
184 <FONT color="green">181</FONT> verifyBracketing(min, max, f);<a name="line.181"></a>
185 <FONT color="green">182</FONT> <a name="line.182"></a>
186 <FONT color="green">183</FONT> double coefficients[] = ((PolynomialFunction) f).getCoefficients();<a name="line.183"></a>
187 <FONT color="green">184</FONT> Complex c[] = new Complex[coefficients.length];<a name="line.184"></a>
188 <FONT color="green">185</FONT> for (int i = 0; i &lt; coefficients.length; i++) {<a name="line.185"></a>
189 <FONT color="green">186</FONT> c[i] = new Complex(coefficients[i], 0.0);<a name="line.186"></a>
190 <FONT color="green">187</FONT> }<a name="line.187"></a>
191 <FONT color="green">188</FONT> Complex initial = new Complex(0.5 * (min + max), 0.0);<a name="line.188"></a>
192 <FONT color="green">189</FONT> Complex z = solve(c, initial);<a name="line.189"></a>
193 <FONT color="green">190</FONT> if (isRootOK(min, max, z)) {<a name="line.190"></a>
194 <FONT color="green">191</FONT> setResult(z.getReal(), iterationCount);<a name="line.191"></a>
195 <FONT color="green">192</FONT> return result;<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> // solve all roots and select the one we're seeking<a name="line.195"></a>
199 <FONT color="green">196</FONT> Complex[] root = solveAll(c, initial);<a name="line.196"></a>
200 <FONT color="green">197</FONT> for (int i = 0; i &lt; root.length; i++) {<a name="line.197"></a>
201 <FONT color="green">198</FONT> if (isRootOK(min, max, root[i])) {<a name="line.198"></a>
202 <FONT color="green">199</FONT> setResult(root[i].getReal(), iterationCount);<a name="line.199"></a>
203 <FONT color="green">200</FONT> return result;<a name="line.200"></a>
204 <FONT color="green">201</FONT> }<a name="line.201"></a>
205 <FONT color="green">202</FONT> }<a name="line.202"></a>
206 <FONT color="green">203</FONT> <a name="line.203"></a>
207 <FONT color="green">204</FONT> // should never happen<a name="line.204"></a>
208 <FONT color="green">205</FONT> throw new ConvergenceException();<a name="line.205"></a>
209 <FONT color="green">206</FONT> }<a name="line.206"></a>
210 <FONT color="green">207</FONT> <a name="line.207"></a>
211 <FONT color="green">208</FONT> /**<a name="line.208"></a>
212 <FONT color="green">209</FONT> * Returns true iff the given complex root is actually a real zero<a name="line.209"></a>
213 <FONT color="green">210</FONT> * in the given interval, within the solver tolerance level.<a name="line.210"></a>
214 <FONT color="green">211</FONT> *<a name="line.211"></a>
215 <FONT color="green">212</FONT> * @param min the lower bound for the interval<a name="line.212"></a>
216 <FONT color="green">213</FONT> * @param max the upper bound for the interval<a name="line.213"></a>
217 <FONT color="green">214</FONT> * @param z the complex root<a name="line.214"></a>
218 <FONT color="green">215</FONT> * @return true iff z is the sought-after real zero<a name="line.215"></a>
219 <FONT color="green">216</FONT> */<a name="line.216"></a>
220 <FONT color="green">217</FONT> protected boolean isRootOK(double min, double max, Complex z) {<a name="line.217"></a>
221 <FONT color="green">218</FONT> double tolerance = Math.max(relativeAccuracy * z.abs(), absoluteAccuracy);<a name="line.218"></a>
222 <FONT color="green">219</FONT> return (isSequence(min, z.getReal(), max)) &amp;&amp;<a name="line.219"></a>
223 <FONT color="green">220</FONT> (Math.abs(z.getImaginary()) &lt;= tolerance ||<a name="line.220"></a>
224 <FONT color="green">221</FONT> z.abs() &lt;= functionValueAccuracy);<a name="line.221"></a>
225 <FONT color="green">222</FONT> }<a name="line.222"></a>
226 <FONT color="green">223</FONT> <a name="line.223"></a>
227 <FONT color="green">224</FONT> /**<a name="line.224"></a>
228 <FONT color="green">225</FONT> * Find all complex roots for the polynomial with the given coefficients,<a name="line.225"></a>
229 <FONT color="green">226</FONT> * starting from the given initial value.<a name="line.226"></a>
230 <FONT color="green">227</FONT> *<a name="line.227"></a>
231 <FONT color="green">228</FONT> * @param coefficients the polynomial coefficients array<a name="line.228"></a>
232 <FONT color="green">229</FONT> * @param initial the start value to use<a name="line.229"></a>
233 <FONT color="green">230</FONT> * @return the point at which the function value is zero<a name="line.230"></a>
234 <FONT color="green">231</FONT> * @throws ConvergenceException if the maximum iteration count is exceeded<a name="line.231"></a>
235 <FONT color="green">232</FONT> * or the solver detects convergence problems otherwise<a name="line.232"></a>
236 <FONT color="green">233</FONT> * @throws FunctionEvaluationException if an error occurs evaluating the<a name="line.233"></a>
237 <FONT color="green">234</FONT> * function<a name="line.234"></a>
238 <FONT color="green">235</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.235"></a>
239 <FONT color="green">236</FONT> */<a name="line.236"></a>
240 <FONT color="green">237</FONT> public Complex[] solveAll(double coefficients[], double initial) throws<a name="line.237"></a>
241 <FONT color="green">238</FONT> ConvergenceException, FunctionEvaluationException {<a name="line.238"></a>
242 <FONT color="green">239</FONT> <a name="line.239"></a>
243 <FONT color="green">240</FONT> Complex c[] = new Complex[coefficients.length];<a name="line.240"></a>
244 <FONT color="green">241</FONT> Complex z = new Complex(initial, 0.0);<a name="line.241"></a>
245 <FONT color="green">242</FONT> for (int i = 0; i &lt; c.length; i++) {<a name="line.242"></a>
246 <FONT color="green">243</FONT> c[i] = new Complex(coefficients[i], 0.0);<a name="line.243"></a>
247 <FONT color="green">244</FONT> }<a name="line.244"></a>
248 <FONT color="green">245</FONT> return solveAll(c, z);<a name="line.245"></a>
249 <FONT color="green">246</FONT> }<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> * Find all complex roots for the polynomial with the given coefficients,<a name="line.249"></a>
253 <FONT color="green">250</FONT> * starting from the given initial value.<a name="line.250"></a>
254 <FONT color="green">251</FONT> *<a name="line.251"></a>
255 <FONT color="green">252</FONT> * @param coefficients the polynomial coefficients array<a name="line.252"></a>
256 <FONT color="green">253</FONT> * @param initial the start value to use<a name="line.253"></a>
257 <FONT color="green">254</FONT> * @return the point at which the function value is zero<a name="line.254"></a>
258 <FONT color="green">255</FONT> * @throws MaxIterationsExceededException if the maximum iteration count is exceeded<a name="line.255"></a>
259 <FONT color="green">256</FONT> * or the solver detects convergence problems otherwise<a name="line.256"></a>
260 <FONT color="green">257</FONT> * @throws FunctionEvaluationException if an error occurs evaluating the<a name="line.257"></a>
261 <FONT color="green">258</FONT> * function<a name="line.258"></a>
262 <FONT color="green">259</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.259"></a>
263 <FONT color="green">260</FONT> */<a name="line.260"></a>
264 <FONT color="green">261</FONT> public Complex[] solveAll(Complex coefficients[], Complex initial) throws<a name="line.261"></a>
265 <FONT color="green">262</FONT> MaxIterationsExceededException, FunctionEvaluationException {<a name="line.262"></a>
266 <FONT color="green">263</FONT> <a name="line.263"></a>
267 <FONT color="green">264</FONT> int n = coefficients.length - 1;<a name="line.264"></a>
268 <FONT color="green">265</FONT> int iterationCount = 0;<a name="line.265"></a>
269 <FONT color="green">266</FONT> if (n &lt; 1) {<a name="line.266"></a>
270 <FONT color="green">267</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.267"></a>
271 <FONT color="green">268</FONT> NON_POSITIVE_DEGREE_MESSAGE, n);<a name="line.268"></a>
272 <FONT color="green">269</FONT> }<a name="line.269"></a>
273 <FONT color="green">270</FONT> Complex c[] = new Complex[n+1]; // coefficients for deflated polynomial<a name="line.270"></a>
274 <FONT color="green">271</FONT> for (int i = 0; i &lt;= n; i++) {<a name="line.271"></a>
275 <FONT color="green">272</FONT> c[i] = coefficients[i];<a name="line.272"></a>
276 <FONT color="green">273</FONT> }<a name="line.273"></a>
277 <FONT color="green">274</FONT> <a name="line.274"></a>
278 <FONT color="green">275</FONT> // solve individual root successively<a name="line.275"></a>
279 <FONT color="green">276</FONT> Complex root[] = new Complex[n];<a name="line.276"></a>
280 <FONT color="green">277</FONT> for (int i = 0; i &lt; n; i++) {<a name="line.277"></a>
281 <FONT color="green">278</FONT> Complex subarray[] = new Complex[n-i+1];<a name="line.278"></a>
282 <FONT color="green">279</FONT> System.arraycopy(c, 0, subarray, 0, subarray.length);<a name="line.279"></a>
283 <FONT color="green">280</FONT> root[i] = solve(subarray, initial);<a name="line.280"></a>
284 <FONT color="green">281</FONT> // polynomial deflation using synthetic division<a name="line.281"></a>
285 <FONT color="green">282</FONT> Complex newc = c[n-i];<a name="line.282"></a>
286 <FONT color="green">283</FONT> Complex oldc = null;<a name="line.283"></a>
287 <FONT color="green">284</FONT> for (int j = n-i-1; j &gt;= 0; j--) {<a name="line.284"></a>
288 <FONT color="green">285</FONT> oldc = c[j];<a name="line.285"></a>
289 <FONT color="green">286</FONT> c[j] = newc;<a name="line.286"></a>
290 <FONT color="green">287</FONT> newc = oldc.add(newc.multiply(root[i]));<a name="line.287"></a>
291 <FONT color="green">288</FONT> }<a name="line.288"></a>
292 <FONT color="green">289</FONT> iterationCount += this.iterationCount;<a name="line.289"></a>
293 <FONT color="green">290</FONT> }<a name="line.290"></a>
294 <FONT color="green">291</FONT> <a name="line.291"></a>
295 <FONT color="green">292</FONT> resultComputed = true;<a name="line.292"></a>
296 <FONT color="green">293</FONT> this.iterationCount = iterationCount;<a name="line.293"></a>
297 <FONT color="green">294</FONT> return root;<a name="line.294"></a>
298 <FONT color="green">295</FONT> }<a name="line.295"></a>
299 <FONT color="green">296</FONT> <a name="line.296"></a>
300 <FONT color="green">297</FONT> /**<a name="line.297"></a>
301 <FONT color="green">298</FONT> * Find a complex root for the polynomial with the given coefficients,<a name="line.298"></a>
302 <FONT color="green">299</FONT> * starting from the given initial value.<a name="line.299"></a>
303 <FONT color="green">300</FONT> *<a name="line.300"></a>
304 <FONT color="green">301</FONT> * @param coefficients the polynomial coefficients array<a name="line.301"></a>
305 <FONT color="green">302</FONT> * @param initial the start value to use<a name="line.302"></a>
306 <FONT color="green">303</FONT> * @return the point at which the function value is zero<a name="line.303"></a>
307 <FONT color="green">304</FONT> * @throws MaxIterationsExceededException if the maximum iteration count is exceeded<a name="line.304"></a>
308 <FONT color="green">305</FONT> * or the solver detects convergence problems otherwise<a name="line.305"></a>
309 <FONT color="green">306</FONT> * @throws FunctionEvaluationException if an error occurs evaluating the<a name="line.306"></a>
310 <FONT color="green">307</FONT> * function<a name="line.307"></a>
311 <FONT color="green">308</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.308"></a>
312 <FONT color="green">309</FONT> */<a name="line.309"></a>
313 <FONT color="green">310</FONT> public Complex solve(Complex coefficients[], Complex initial) throws<a name="line.310"></a>
314 <FONT color="green">311</FONT> MaxIterationsExceededException, FunctionEvaluationException {<a name="line.311"></a>
315 <FONT color="green">312</FONT> <a name="line.312"></a>
316 <FONT color="green">313</FONT> int n = coefficients.length - 1;<a name="line.313"></a>
317 <FONT color="green">314</FONT> if (n &lt; 1) {<a name="line.314"></a>
318 <FONT color="green">315</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.315"></a>
319 <FONT color="green">316</FONT> NON_POSITIVE_DEGREE_MESSAGE, n);<a name="line.316"></a>
320 <FONT color="green">317</FONT> }<a name="line.317"></a>
321 <FONT color="green">318</FONT> Complex N = new Complex(n, 0.0);<a name="line.318"></a>
322 <FONT color="green">319</FONT> Complex N1 = new Complex(n - 1, 0.0);<a name="line.319"></a>
323 <FONT color="green">320</FONT> <a name="line.320"></a>
324 <FONT color="green">321</FONT> int i = 1;<a name="line.321"></a>
325 <FONT color="green">322</FONT> Complex pv = null;<a name="line.322"></a>
326 <FONT color="green">323</FONT> Complex dv = null;<a name="line.323"></a>
327 <FONT color="green">324</FONT> Complex d2v = null;<a name="line.324"></a>
328 <FONT color="green">325</FONT> Complex G = null;<a name="line.325"></a>
329 <FONT color="green">326</FONT> Complex G2 = null;<a name="line.326"></a>
330 <FONT color="green">327</FONT> Complex H = null;<a name="line.327"></a>
331 <FONT color="green">328</FONT> Complex delta = null;<a name="line.328"></a>
332 <FONT color="green">329</FONT> Complex denominator = null;<a name="line.329"></a>
333 <FONT color="green">330</FONT> Complex z = initial;<a name="line.330"></a>
334 <FONT color="green">331</FONT> Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);<a name="line.331"></a>
335 <FONT color="green">332</FONT> while (i &lt;= maximalIterationCount) {<a name="line.332"></a>
336 <FONT color="green">333</FONT> // Compute pv (polynomial value), dv (derivative value), and<a name="line.333"></a>
337 <FONT color="green">334</FONT> // d2v (second derivative value) simultaneously.<a name="line.334"></a>
338 <FONT color="green">335</FONT> pv = coefficients[n];<a name="line.335"></a>
339 <FONT color="green">336</FONT> dv = Complex.ZERO;<a name="line.336"></a>
340 <FONT color="green">337</FONT> d2v = Complex.ZERO;<a name="line.337"></a>
341 <FONT color="green">338</FONT> for (int j = n-1; j &gt;= 0; j--) {<a name="line.338"></a>
342 <FONT color="green">339</FONT> d2v = dv.add(z.multiply(d2v));<a name="line.339"></a>
343 <FONT color="green">340</FONT> dv = pv.add(z.multiply(dv));<a name="line.340"></a>
344 <FONT color="green">341</FONT> pv = coefficients[j].add(z.multiply(pv));<a name="line.341"></a>
345 <FONT color="green">342</FONT> }<a name="line.342"></a>
346 <FONT color="green">343</FONT> d2v = d2v.multiply(new Complex(2.0, 0.0));<a name="line.343"></a>
347 <FONT color="green">344</FONT> <a name="line.344"></a>
348 <FONT color="green">345</FONT> // check for convergence<a name="line.345"></a>
349 <FONT color="green">346</FONT> double tolerance = Math.max(relativeAccuracy * z.abs(),<a name="line.346"></a>
350 <FONT color="green">347</FONT> absoluteAccuracy);<a name="line.347"></a>
351 <FONT color="green">348</FONT> if ((z.subtract(oldz)).abs() &lt;= tolerance) {<a name="line.348"></a>
352 <FONT color="green">349</FONT> resultComputed = true;<a name="line.349"></a>
353 <FONT color="green">350</FONT> iterationCount = i;<a name="line.350"></a>
354 <FONT color="green">351</FONT> return z;<a name="line.351"></a>
355 <FONT color="green">352</FONT> }<a name="line.352"></a>
356 <FONT color="green">353</FONT> if (pv.abs() &lt;= functionValueAccuracy) {<a name="line.353"></a>
357 <FONT color="green">354</FONT> resultComputed = true;<a name="line.354"></a>
358 <FONT color="green">355</FONT> iterationCount = i;<a name="line.355"></a>
359 <FONT color="green">356</FONT> return z;<a name="line.356"></a>
360 <FONT color="green">357</FONT> }<a name="line.357"></a>
361 <FONT color="green">358</FONT> <a name="line.358"></a>
362 <FONT color="green">359</FONT> // now pv != 0, calculate the new approximation<a name="line.359"></a>
363 <FONT color="green">360</FONT> G = dv.divide(pv);<a name="line.360"></a>
364 <FONT color="green">361</FONT> G2 = G.multiply(G);<a name="line.361"></a>
365 <FONT color="green">362</FONT> H = G2.subtract(d2v.divide(pv));<a name="line.362"></a>
366 <FONT color="green">363</FONT> delta = N1.multiply((N.multiply(H)).subtract(G2));<a name="line.363"></a>
367 <FONT color="green">364</FONT> // choose a denominator larger in magnitude<a name="line.364"></a>
368 <FONT color="green">365</FONT> Complex deltaSqrt = delta.sqrt();<a name="line.365"></a>
369 <FONT color="green">366</FONT> Complex dplus = G.add(deltaSqrt);<a name="line.366"></a>
370 <FONT color="green">367</FONT> Complex dminus = G.subtract(deltaSqrt);<a name="line.367"></a>
371 <FONT color="green">368</FONT> denominator = dplus.abs() &gt; dminus.abs() ? dplus : dminus;<a name="line.368"></a>
372 <FONT color="green">369</FONT> // Perturb z if denominator is zero, for instance,<a name="line.369"></a>
373 <FONT color="green">370</FONT> // p(x) = x^3 + 1, z = 0.<a name="line.370"></a>
374 <FONT color="green">371</FONT> if (denominator.equals(new Complex(0.0, 0.0))) {<a name="line.371"></a>
375 <FONT color="green">372</FONT> z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));<a name="line.372"></a>
376 <FONT color="green">373</FONT> oldz = new Complex(Double.POSITIVE_INFINITY,<a name="line.373"></a>
377 <FONT color="green">374</FONT> Double.POSITIVE_INFINITY);<a name="line.374"></a>
378 <FONT color="green">375</FONT> } else {<a name="line.375"></a>
379 <FONT color="green">376</FONT> oldz = z;<a name="line.376"></a>
380 <FONT color="green">377</FONT> z = z.subtract(N.divide(denominator));<a name="line.377"></a>
381 <FONT color="green">378</FONT> }<a name="line.378"></a>
382 <FONT color="green">379</FONT> i++;<a name="line.379"></a>
383 <FONT color="green">380</FONT> }<a name="line.380"></a>
384 <FONT color="green">381</FONT> throw new MaxIterationsExceededException(maximalIterationCount);<a name="line.381"></a>
385 <FONT color="green">382</FONT> }<a name="line.382"></a>
386 <FONT color="green">383</FONT> }<a name="line.383"></a>
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447 </PRE>
448 </BODY>
449 </HTML>