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

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.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.polynomials;<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.DuplicateSampleAbscissaException;<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.analysis.UnivariateRealFunction;<a name="line.22"></a>
26 <FONT color="green">023</FONT> <a name="line.23"></a>
27 <FONT color="green">024</FONT> /**<a name="line.24"></a>
28 <FONT color="green">025</FONT> * Implements the representation of a real polynomial function in<a name="line.25"></a>
29 <FONT color="green">026</FONT> * &lt;a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html"&gt;<a name="line.26"></a>
30 <FONT color="green">027</FONT> * Lagrange Form&lt;/a&gt;. For reference, see &lt;b&gt;Introduction to Numerical<a name="line.27"></a>
31 <FONT color="green">028</FONT> * Analysis&lt;/b&gt;, ISBN 038795452X, chapter 2.<a name="line.28"></a>
32 <FONT color="green">029</FONT> * &lt;p&gt;<a name="line.29"></a>
33 <FONT color="green">030</FONT> * The approximated function should be smooth enough for Lagrange polynomial<a name="line.30"></a>
34 <FONT color="green">031</FONT> * to work well. Otherwise, consider using splines instead.&lt;/p&gt;<a name="line.31"></a>
35 <FONT color="green">032</FONT> *<a name="line.32"></a>
36 <FONT color="green">033</FONT> * @version $Revision: 922708 $ $Date: 2010-03-13 20:15:47 -0500 (Sat, 13 Mar 2010) $<a name="line.33"></a>
37 <FONT color="green">034</FONT> * @since 1.2<a name="line.34"></a>
38 <FONT color="green">035</FONT> */<a name="line.35"></a>
39 <FONT color="green">036</FONT> public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {<a name="line.36"></a>
40 <FONT color="green">037</FONT> <a name="line.37"></a>
41 <FONT color="green">038</FONT> /**<a name="line.38"></a>
42 <FONT color="green">039</FONT> * The coefficients of the polynomial, ordered by degree -- i.e.<a name="line.39"></a>
43 <FONT color="green">040</FONT> * coefficients[0] is the constant term and coefficients[n] is the<a name="line.40"></a>
44 <FONT color="green">041</FONT> * coefficient of x^n where n is the degree of the polynomial.<a name="line.41"></a>
45 <FONT color="green">042</FONT> */<a name="line.42"></a>
46 <FONT color="green">043</FONT> private double coefficients[];<a name="line.43"></a>
47 <FONT color="green">044</FONT> <a name="line.44"></a>
48 <FONT color="green">045</FONT> /**<a name="line.45"></a>
49 <FONT color="green">046</FONT> * Interpolating points (abscissas).<a name="line.46"></a>
50 <FONT color="green">047</FONT> */<a name="line.47"></a>
51 <FONT color="green">048</FONT> private final double x[];<a name="line.48"></a>
52 <FONT color="green">049</FONT> <a name="line.49"></a>
53 <FONT color="green">050</FONT> /**<a name="line.50"></a>
54 <FONT color="green">051</FONT> * Function values at interpolating points.<a name="line.51"></a>
55 <FONT color="green">052</FONT> */<a name="line.52"></a>
56 <FONT color="green">053</FONT> private final double y[];<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> * Whether the polynomial coefficients are available.<a name="line.56"></a>
60 <FONT color="green">057</FONT> */<a name="line.57"></a>
61 <FONT color="green">058</FONT> private boolean coefficientsComputed;<a name="line.58"></a>
62 <FONT color="green">059</FONT> <a name="line.59"></a>
63 <FONT color="green">060</FONT> /**<a name="line.60"></a>
64 <FONT color="green">061</FONT> * Construct a Lagrange polynomial with the given abscissas and function<a name="line.61"></a>
65 <FONT color="green">062</FONT> * values. The order of interpolating points are not important.<a name="line.62"></a>
66 <FONT color="green">063</FONT> * &lt;p&gt;<a name="line.63"></a>
67 <FONT color="green">064</FONT> * The constructor makes copy of the input arrays and assigns them.&lt;/p&gt;<a name="line.64"></a>
68 <FONT color="green">065</FONT> *<a name="line.65"></a>
69 <FONT color="green">066</FONT> * @param x interpolating points<a name="line.66"></a>
70 <FONT color="green">067</FONT> * @param y function values at interpolating points<a name="line.67"></a>
71 <FONT color="green">068</FONT> * @throws IllegalArgumentException if input arrays are not valid<a name="line.68"></a>
72 <FONT color="green">069</FONT> */<a name="line.69"></a>
73 <FONT color="green">070</FONT> public PolynomialFunctionLagrangeForm(double x[], double y[])<a name="line.70"></a>
74 <FONT color="green">071</FONT> throws IllegalArgumentException {<a name="line.71"></a>
75 <FONT color="green">072</FONT> <a name="line.72"></a>
76 <FONT color="green">073</FONT> verifyInterpolationArray(x, y);<a name="line.73"></a>
77 <FONT color="green">074</FONT> this.x = new double[x.length];<a name="line.74"></a>
78 <FONT color="green">075</FONT> this.y = new double[y.length];<a name="line.75"></a>
79 <FONT color="green">076</FONT> System.arraycopy(x, 0, this.x, 0, x.length);<a name="line.76"></a>
80 <FONT color="green">077</FONT> System.arraycopy(y, 0, this.y, 0, y.length);<a name="line.77"></a>
81 <FONT color="green">078</FONT> coefficientsComputed = false;<a name="line.78"></a>
82 <FONT color="green">079</FONT> }<a name="line.79"></a>
83 <FONT color="green">080</FONT> <a name="line.80"></a>
84 <FONT color="green">081</FONT> /**<a name="line.81"></a>
85 <FONT color="green">082</FONT> * Calculate the function value at the given point.<a name="line.82"></a>
86 <FONT color="green">083</FONT> *<a name="line.83"></a>
87 <FONT color="green">084</FONT> * @param z the point at which the function value is to be computed<a name="line.84"></a>
88 <FONT color="green">085</FONT> * @return the function value<a name="line.85"></a>
89 <FONT color="green">086</FONT> * @throws FunctionEvaluationException if a runtime error occurs<a name="line.86"></a>
90 <FONT color="green">087</FONT> * @see UnivariateRealFunction#value(double)<a name="line.87"></a>
91 <FONT color="green">088</FONT> */<a name="line.88"></a>
92 <FONT color="green">089</FONT> public double value(double z) throws FunctionEvaluationException {<a name="line.89"></a>
93 <FONT color="green">090</FONT> try {<a name="line.90"></a>
94 <FONT color="green">091</FONT> return evaluate(x, y, z);<a name="line.91"></a>
95 <FONT color="green">092</FONT> } catch (DuplicateSampleAbscissaException e) {<a name="line.92"></a>
96 <FONT color="green">093</FONT> throw new FunctionEvaluationException(e, z, e.getPattern(), e.getArguments());<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> <a name="line.96"></a>
100 <FONT color="green">097</FONT> /**<a name="line.97"></a>
101 <FONT color="green">098</FONT> * Returns the degree of the polynomial.<a name="line.98"></a>
102 <FONT color="green">099</FONT> *<a name="line.99"></a>
103 <FONT color="green">100</FONT> * @return the degree of the polynomial<a name="line.100"></a>
104 <FONT color="green">101</FONT> */<a name="line.101"></a>
105 <FONT color="green">102</FONT> public int degree() {<a name="line.102"></a>
106 <FONT color="green">103</FONT> return x.length - 1;<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> * Returns a copy of the interpolating points array.<a name="line.107"></a>
111 <FONT color="green">108</FONT> * &lt;p&gt;<a name="line.108"></a>
112 <FONT color="green">109</FONT> * Changes made to the returned copy will not affect the polynomial.&lt;/p&gt;<a name="line.109"></a>
113 <FONT color="green">110</FONT> *<a name="line.110"></a>
114 <FONT color="green">111</FONT> * @return a fresh copy of the interpolating points array<a name="line.111"></a>
115 <FONT color="green">112</FONT> */<a name="line.112"></a>
116 <FONT color="green">113</FONT> public double[] getInterpolatingPoints() {<a name="line.113"></a>
117 <FONT color="green">114</FONT> double[] out = new double[x.length];<a name="line.114"></a>
118 <FONT color="green">115</FONT> System.arraycopy(x, 0, out, 0, x.length);<a name="line.115"></a>
119 <FONT color="green">116</FONT> return out;<a name="line.116"></a>
120 <FONT color="green">117</FONT> }<a name="line.117"></a>
121 <FONT color="green">118</FONT> <a name="line.118"></a>
122 <FONT color="green">119</FONT> /**<a name="line.119"></a>
123 <FONT color="green">120</FONT> * Returns a copy of the interpolating values array.<a name="line.120"></a>
124 <FONT color="green">121</FONT> * &lt;p&gt;<a name="line.121"></a>
125 <FONT color="green">122</FONT> * Changes made to the returned copy will not affect the polynomial.&lt;/p&gt;<a name="line.122"></a>
126 <FONT color="green">123</FONT> *<a name="line.123"></a>
127 <FONT color="green">124</FONT> * @return a fresh copy of the interpolating values array<a name="line.124"></a>
128 <FONT color="green">125</FONT> */<a name="line.125"></a>
129 <FONT color="green">126</FONT> public double[] getInterpolatingValues() {<a name="line.126"></a>
130 <FONT color="green">127</FONT> double[] out = new double[y.length];<a name="line.127"></a>
131 <FONT color="green">128</FONT> System.arraycopy(y, 0, out, 0, y.length);<a name="line.128"></a>
132 <FONT color="green">129</FONT> return out;<a name="line.129"></a>
133 <FONT color="green">130</FONT> }<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> * Returns a copy of the coefficients array.<a name="line.133"></a>
137 <FONT color="green">134</FONT> * &lt;p&gt;<a name="line.134"></a>
138 <FONT color="green">135</FONT> * Changes made to the returned copy will not affect the polynomial.&lt;/p&gt;<a name="line.135"></a>
139 <FONT color="green">136</FONT> * &lt;p&gt;<a name="line.136"></a>
140 <FONT color="green">137</FONT> * Note that coefficients computation can be ill-conditioned. Use with caution<a name="line.137"></a>
141 <FONT color="green">138</FONT> * and only when it is necessary.&lt;/p&gt;<a name="line.138"></a>
142 <FONT color="green">139</FONT> *<a name="line.139"></a>
143 <FONT color="green">140</FONT> * @return a fresh copy of the coefficients array<a name="line.140"></a>
144 <FONT color="green">141</FONT> */<a name="line.141"></a>
145 <FONT color="green">142</FONT> public double[] getCoefficients() {<a name="line.142"></a>
146 <FONT color="green">143</FONT> if (!coefficientsComputed) {<a name="line.143"></a>
147 <FONT color="green">144</FONT> computeCoefficients();<a name="line.144"></a>
148 <FONT color="green">145</FONT> }<a name="line.145"></a>
149 <FONT color="green">146</FONT> double[] out = new double[coefficients.length];<a name="line.146"></a>
150 <FONT color="green">147</FONT> System.arraycopy(coefficients, 0, out, 0, coefficients.length);<a name="line.147"></a>
151 <FONT color="green">148</FONT> return out;<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> /**<a name="line.151"></a>
155 <FONT color="green">152</FONT> * Evaluate the Lagrange polynomial using<a name="line.152"></a>
156 <FONT color="green">153</FONT> * &lt;a href="http://mathworld.wolfram.com/NevillesAlgorithm.html"&gt;<a name="line.153"></a>
157 <FONT color="green">154</FONT> * Neville's Algorithm&lt;/a&gt;. It takes O(N^2) time.<a name="line.154"></a>
158 <FONT color="green">155</FONT> * &lt;p&gt;<a name="line.155"></a>
159 <FONT color="green">156</FONT> * This function is made public static so that users can call it directly<a name="line.156"></a>
160 <FONT color="green">157</FONT> * without instantiating PolynomialFunctionLagrangeForm object.&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 x the interpolating points array<a name="line.159"></a>
163 <FONT color="green">160</FONT> * @param y the interpolating values array<a name="line.160"></a>
164 <FONT color="green">161</FONT> * @param z the point at which the function value is to be computed<a name="line.161"></a>
165 <FONT color="green">162</FONT> * @return the function value<a name="line.162"></a>
166 <FONT color="green">163</FONT> * @throws DuplicateSampleAbscissaException if the sample has duplicate abscissas<a name="line.163"></a>
167 <FONT color="green">164</FONT> * @throws IllegalArgumentException if inputs are not valid<a name="line.164"></a>
168 <FONT color="green">165</FONT> */<a name="line.165"></a>
169 <FONT color="green">166</FONT> public static double evaluate(double x[], double y[], double z) throws<a name="line.166"></a>
170 <FONT color="green">167</FONT> DuplicateSampleAbscissaException, IllegalArgumentException {<a name="line.167"></a>
171 <FONT color="green">168</FONT> <a name="line.168"></a>
172 <FONT color="green">169</FONT> verifyInterpolationArray(x, y);<a name="line.169"></a>
173 <FONT color="green">170</FONT> <a name="line.170"></a>
174 <FONT color="green">171</FONT> int nearest = 0;<a name="line.171"></a>
175 <FONT color="green">172</FONT> final int n = x.length;<a name="line.172"></a>
176 <FONT color="green">173</FONT> final double[] c = new double[n];<a name="line.173"></a>
177 <FONT color="green">174</FONT> final double[] d = new double[n];<a name="line.174"></a>
178 <FONT color="green">175</FONT> double min_dist = Double.POSITIVE_INFINITY;<a name="line.175"></a>
179 <FONT color="green">176</FONT> for (int i = 0; i &lt; n; i++) {<a name="line.176"></a>
180 <FONT color="green">177</FONT> // initialize the difference arrays<a name="line.177"></a>
181 <FONT color="green">178</FONT> c[i] = y[i];<a name="line.178"></a>
182 <FONT color="green">179</FONT> d[i] = y[i];<a name="line.179"></a>
183 <FONT color="green">180</FONT> // find out the abscissa closest to z<a name="line.180"></a>
184 <FONT color="green">181</FONT> final double dist = Math.abs(z - x[i]);<a name="line.181"></a>
185 <FONT color="green">182</FONT> if (dist &lt; min_dist) {<a name="line.182"></a>
186 <FONT color="green">183</FONT> nearest = i;<a name="line.183"></a>
187 <FONT color="green">184</FONT> min_dist = dist;<a name="line.184"></a>
188 <FONT color="green">185</FONT> }<a name="line.185"></a>
189 <FONT color="green">186</FONT> }<a name="line.186"></a>
190 <FONT color="green">187</FONT> <a name="line.187"></a>
191 <FONT color="green">188</FONT> // initial approximation to the function value at z<a name="line.188"></a>
192 <FONT color="green">189</FONT> double value = y[nearest];<a name="line.189"></a>
193 <FONT color="green">190</FONT> <a name="line.190"></a>
194 <FONT color="green">191</FONT> for (int i = 1; i &lt; n; i++) {<a name="line.191"></a>
195 <FONT color="green">192</FONT> for (int j = 0; j &lt; n-i; j++) {<a name="line.192"></a>
196 <FONT color="green">193</FONT> final double tc = x[j] - z;<a name="line.193"></a>
197 <FONT color="green">194</FONT> final double td = x[i+j] - z;<a name="line.194"></a>
198 <FONT color="green">195</FONT> final double divider = x[j] - x[i+j];<a name="line.195"></a>
199 <FONT color="green">196</FONT> if (divider == 0.0) {<a name="line.196"></a>
200 <FONT color="green">197</FONT> // This happens only when two abscissas are identical.<a name="line.197"></a>
201 <FONT color="green">198</FONT> throw new DuplicateSampleAbscissaException(x[i], i, i+j);<a name="line.198"></a>
202 <FONT color="green">199</FONT> }<a name="line.199"></a>
203 <FONT color="green">200</FONT> // update the difference arrays<a name="line.200"></a>
204 <FONT color="green">201</FONT> final double w = (c[j+1] - d[j]) / divider;<a name="line.201"></a>
205 <FONT color="green">202</FONT> c[j] = tc * w;<a name="line.202"></a>
206 <FONT color="green">203</FONT> d[j] = td * w;<a name="line.203"></a>
207 <FONT color="green">204</FONT> }<a name="line.204"></a>
208 <FONT color="green">205</FONT> // sum up the difference terms to get the final value<a name="line.205"></a>
209 <FONT color="green">206</FONT> if (nearest &lt; 0.5*(n-i+1)) {<a name="line.206"></a>
210 <FONT color="green">207</FONT> value += c[nearest]; // fork down<a name="line.207"></a>
211 <FONT color="green">208</FONT> } else {<a name="line.208"></a>
212 <FONT color="green">209</FONT> nearest--;<a name="line.209"></a>
213 <FONT color="green">210</FONT> value += d[nearest]; // fork up<a name="line.210"></a>
214 <FONT color="green">211</FONT> }<a name="line.211"></a>
215 <FONT color="green">212</FONT> }<a name="line.212"></a>
216 <FONT color="green">213</FONT> <a name="line.213"></a>
217 <FONT color="green">214</FONT> return value;<a name="line.214"></a>
218 <FONT color="green">215</FONT> }<a name="line.215"></a>
219 <FONT color="green">216</FONT> <a name="line.216"></a>
220 <FONT color="green">217</FONT> /**<a name="line.217"></a>
221 <FONT color="green">218</FONT> * Calculate the coefficients of Lagrange polynomial from the<a name="line.218"></a>
222 <FONT color="green">219</FONT> * interpolation data. It takes O(N^2) time.<a name="line.219"></a>
223 <FONT color="green">220</FONT> * &lt;p&gt;<a name="line.220"></a>
224 <FONT color="green">221</FONT> * Note this computation can be ill-conditioned. Use with caution<a name="line.221"></a>
225 <FONT color="green">222</FONT> * and only when it is necessary.&lt;/p&gt;<a name="line.222"></a>
226 <FONT color="green">223</FONT> *<a name="line.223"></a>
227 <FONT color="green">224</FONT> * @throws ArithmeticException if any abscissas coincide<a name="line.224"></a>
228 <FONT color="green">225</FONT> */<a name="line.225"></a>
229 <FONT color="green">226</FONT> protected void computeCoefficients() throws ArithmeticException {<a name="line.226"></a>
230 <FONT color="green">227</FONT> <a name="line.227"></a>
231 <FONT color="green">228</FONT> final int n = degree() + 1;<a name="line.228"></a>
232 <FONT color="green">229</FONT> coefficients = new double[n];<a name="line.229"></a>
233 <FONT color="green">230</FONT> for (int i = 0; i &lt; n; i++) {<a name="line.230"></a>
234 <FONT color="green">231</FONT> coefficients[i] = 0.0;<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> // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])<a name="line.234"></a>
238 <FONT color="green">235</FONT> final double[] c = new double[n+1];<a name="line.235"></a>
239 <FONT color="green">236</FONT> c[0] = 1.0;<a name="line.236"></a>
240 <FONT color="green">237</FONT> for (int i = 0; i &lt; n; i++) {<a name="line.237"></a>
241 <FONT color="green">238</FONT> for (int j = i; j &gt; 0; j--) {<a name="line.238"></a>
242 <FONT color="green">239</FONT> c[j] = c[j-1] - c[j] * x[i];<a name="line.239"></a>
243 <FONT color="green">240</FONT> }<a name="line.240"></a>
244 <FONT color="green">241</FONT> c[0] *= -x[i];<a name="line.241"></a>
245 <FONT color="green">242</FONT> c[i+1] = 1;<a name="line.242"></a>
246 <FONT color="green">243</FONT> }<a name="line.243"></a>
247 <FONT color="green">244</FONT> <a name="line.244"></a>
248 <FONT color="green">245</FONT> final double[] tc = new double[n];<a name="line.245"></a>
249 <FONT color="green">246</FONT> for (int i = 0; i &lt; n; i++) {<a name="line.246"></a>
250 <FONT color="green">247</FONT> // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])<a name="line.247"></a>
251 <FONT color="green">248</FONT> double d = 1;<a name="line.248"></a>
252 <FONT color="green">249</FONT> for (int j = 0; j &lt; n; j++) {<a name="line.249"></a>
253 <FONT color="green">250</FONT> if (i != j) {<a name="line.250"></a>
254 <FONT color="green">251</FONT> d *= x[i] - x[j];<a name="line.251"></a>
255 <FONT color="green">252</FONT> }<a name="line.252"></a>
256 <FONT color="green">253</FONT> }<a name="line.253"></a>
257 <FONT color="green">254</FONT> if (d == 0.0) {<a name="line.254"></a>
258 <FONT color="green">255</FONT> // This happens only when two abscissas are identical.<a name="line.255"></a>
259 <FONT color="green">256</FONT> for (int k = 0; k &lt; n; ++k) {<a name="line.256"></a>
260 <FONT color="green">257</FONT> if ((i != k) &amp;&amp; (x[i] == x[k])) {<a name="line.257"></a>
261 <FONT color="green">258</FONT> throw MathRuntimeException.createArithmeticException("identical abscissas x[{0}] == x[{1}] == {2} cause division by zero",<a name="line.258"></a>
262 <FONT color="green">259</FONT> i, k, x[i]);<a name="line.259"></a>
263 <FONT color="green">260</FONT> }<a name="line.260"></a>
264 <FONT color="green">261</FONT> }<a name="line.261"></a>
265 <FONT color="green">262</FONT> }<a name="line.262"></a>
266 <FONT color="green">263</FONT> final double t = y[i] / d;<a name="line.263"></a>
267 <FONT color="green">264</FONT> // Lagrange polynomial is the sum of n terms, each of which is a<a name="line.264"></a>
268 <FONT color="green">265</FONT> // polynomial of degree n-1. tc[] are the coefficients of the i-th<a name="line.265"></a>
269 <FONT color="green">266</FONT> // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).<a name="line.266"></a>
270 <FONT color="green">267</FONT> tc[n-1] = c[n]; // actually c[n] = 1<a name="line.267"></a>
271 <FONT color="green">268</FONT> coefficients[n-1] += t * tc[n-1];<a name="line.268"></a>
272 <FONT color="green">269</FONT> for (int j = n-2; j &gt;= 0; j--) {<a name="line.269"></a>
273 <FONT color="green">270</FONT> tc[j] = c[j+1] + tc[j+1] * x[i];<a name="line.270"></a>
274 <FONT color="green">271</FONT> coefficients[j] += t * tc[j];<a name="line.271"></a>
275 <FONT color="green">272</FONT> }<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> coefficientsComputed = true;<a name="line.275"></a>
279 <FONT color="green">276</FONT> }<a name="line.276"></a>
280 <FONT color="green">277</FONT> <a name="line.277"></a>
281 <FONT color="green">278</FONT> /**<a name="line.278"></a>
282 <FONT color="green">279</FONT> * Verifies that the interpolation arrays are valid.<a name="line.279"></a>
283 <FONT color="green">280</FONT> * &lt;p&gt;<a name="line.280"></a>
284 <FONT color="green">281</FONT> * The arrays features checked by this method are that both arrays have the<a name="line.281"></a>
285 <FONT color="green">282</FONT> * same length and this length is at least 2.<a name="line.282"></a>
286 <FONT color="green">283</FONT> * &lt;/p&gt;<a name="line.283"></a>
287 <FONT color="green">284</FONT> * &lt;p&gt;<a name="line.284"></a>
288 <FONT color="green">285</FONT> * The interpolating points must be distinct. However it is not<a name="line.285"></a>
289 <FONT color="green">286</FONT> * verified here, it is checked in evaluate() and computeCoefficients().<a name="line.286"></a>
290 <FONT color="green">287</FONT> * &lt;/p&gt;<a name="line.287"></a>
291 <FONT color="green">288</FONT> *<a name="line.288"></a>
292 <FONT color="green">289</FONT> * @param x the interpolating points array<a name="line.289"></a>
293 <FONT color="green">290</FONT> * @param y the interpolating values array<a name="line.290"></a>
294 <FONT color="green">291</FONT> * @throws IllegalArgumentException if not valid<a name="line.291"></a>
295 <FONT color="green">292</FONT> * @see #evaluate(double[], double[], double)<a name="line.292"></a>
296 <FONT color="green">293</FONT> * @see #computeCoefficients()<a name="line.293"></a>
297 <FONT color="green">294</FONT> */<a name="line.294"></a>
298 <FONT color="green">295</FONT> public static void verifyInterpolationArray(double x[], double y[])<a name="line.295"></a>
299 <FONT color="green">296</FONT> throws IllegalArgumentException {<a name="line.296"></a>
300 <FONT color="green">297</FONT> <a name="line.297"></a>
301 <FONT color="green">298</FONT> if (x.length != y.length) {<a name="line.298"></a>
302 <FONT color="green">299</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.299"></a>
303 <FONT color="green">300</FONT> "dimension mismatch {0} != {1}", x.length, y.length);<a name="line.300"></a>
304 <FONT color="green">301</FONT> }<a name="line.301"></a>
305 <FONT color="green">302</FONT> <a name="line.302"></a>
306 <FONT color="green">303</FONT> if (x.length &lt; 2) {<a name="line.303"></a>
307 <FONT color="green">304</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.304"></a>
308 <FONT color="green">305</FONT> "{0} points are required, got only {1}", 2, x.length);<a name="line.305"></a>
309 <FONT color="green">306</FONT> }<a name="line.306"></a>
310 <FONT color="green">307</FONT> <a name="line.307"></a>
311 <FONT color="green">308</FONT> }<a name="line.308"></a>
312 <FONT color="green">309</FONT> }<a name="line.309"></a>
313
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
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373 </PRE>
374 </BODY>
375 </HTML>