# Calculating lagrange polynomial for 100 points?

I need to calculate the lagrange polynomial which approximates $e^x$ at $101$ points, the points $\frac{k}{101^2}$ for $k\in\{0,1,2\dots 100\}$.

I tried the following code:

import java.math.*;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunctionLagrangeForm;

public class interpolation {

public static void main(String[] args) {
// TODO Auto-generated method stub
int N = 100;

double [] X = new double [N+1];
double [] Y = new double [N+1];

// Here we give the computer the points:

for(int k=0; k<=N ;k++){
double a = k/(N*N*1.0);
double b = Math.exp(a);
X[k]=a;
Y[k]=b;
}

// here we build the polynomial

PolynomialFunctionLagrangeForm pol = new PolynomialFunctionLagrangeForm( X , Y );

double [] P = pol.getCoefficients();

// here we print the polynomial

for(int i=P.length -1 ;i >= 0;i--){
System.out.print( P[i] + "    ");
}

}

}


It gives good results for $N=5,10$ but bad ones for $N=50,100$ (it doesn't resemble the Taylor expansion at all). I would like help in calculating this polynomial more precisely. Thank you kindly.

• Are you aware of Runge's phenomenon? High-order polynomial interpolation at evenly spaced points is famously unstable; you'd get much better results using Chebyshev points instead. – Christian Clason Mar 13 '16 at 20:33
• If you are required to solve this for $N=50,100$ for school or something, you are likely expected to see it fail. If you could choose some other interpolation, I would recommend something like Akima splines. – spektr Mar 13 '16 at 23:57
• Actually the optimal approximation of the exponential for small arguments using polynomials / rationals has a closed form solution due to Padé. See en.wikipedia.org/wiki/Pad%C3%A9_approximant – Richard Zhang Mar 14 '16 at 0:30
• @nicoguaro That's strange because in my tests $e^x$ on $[0,10]$ only needs about 30 Chebyshev points to interpolate with almost full double precision, quite easily. The last 70 should be almost zero, on order of $100e^{10}\epsilon$. – Kirill Mar 14 '16 at 22:47
• @nicoguaro I don't see it, it works just fine with Chebyshev points: x = (cos((2*linspace(1,n,n)-1)*pi/(2*n))+1)*5; fb = BarycentricInterpolator(x, exp(x)) (interp1d fails very badly with rel. errors like $0.003$, otoh). Indeed, it should work fine, based on the mathematics; polynomial interpolation with Chebyshev points is equivalent to the discrete cosine transform of the data, see also Berrut-Trefethen's article on polynomial interpolation. Also, you can see that chebfun works just by using it. – Kirill Mar 14 '16 at 23:12