Lagrange Interpolation: Investigating Interpolation Points with C++/Fortran

  • Context: C/C++ 
  • Thread starter Thread starter laminatedevildoll
  • Start date Start date
  • Tags Tags
    Interpolation Lagrange
Click For Summary

Discussion Overview

The discussion revolves around implementing n-point Lagrange interpolation using C++ and Fortran, specifically investigating the quality of interpolation points for the function f(x)=sin(x^2). Participants explore how to write the program, the selection of interpolation points, and the mathematical formulation of Lagrange coefficients.

Discussion Character

  • Technical explanation
  • Homework-related
  • Mathematical reasoning
  • Debate/contested

Main Points Raised

  • One participant expresses uncertainty about how to start writing a program for Lagrange interpolation given the function f(x)=sin(x^2) and the number of points n.
  • Another participant explains the process of calculating Lagrange coefficients and suggests a basic algorithm for implementing the interpolation.
  • There is a discussion about whether the user can select the value of n, with some participants confirming that n can be adjusted by the user.
  • Participants discuss the calculation of x-coordinates based on the value of n and the endpoints of the interval, suggesting a function to simplify this process.
  • A participant shares their code but encounters issues, prompting others to point out potential mistakes, including the incorrect function definition in the code.
  • Suggestions are made to evaluate the polynomial at multiple points to investigate the quality of the interpolation and to experiment with different interpolation points.

Areas of Agreement / Disagreement

Participants generally agree on the process of implementing Lagrange interpolation and the importance of correctly defining the function. However, there are unresolved issues regarding the implementation details in the code and the selection of interpolation points.

Contextual Notes

Some limitations in the discussion include unclear assumptions about the input parameters and the handling of edge cases in the code. There are also unresolved mathematical steps related to the implementation of the algorithm.

laminatedevildoll
Messages
211
Reaction score
0
In my comp physics class, we've been introduced to both c++ and fortran languages. For instance, for our first assignment, I am not sure how to go about investigating the quality of interpolation points for i.e f(x)=sin(x^2) by using n-point langrange interpolation, where n is an input parameter. On paper, I've already calculated P(x) (lagrange) for 10 points in the interval [0,5]. But, how should I start writing the program?
 
Technology news on Phys.org
Are you given the points or are you allowed to obtain the points with some strategy? If you have n points, then computing the Lagrange coefficients and interpolating polynomial, f(x), is not very hard. Mostly you should just follow the mathematical expression for the Lagrange coefficients. From what i remember if you have 3 points (x1,y1), (x2,y2), (x3,y3), then you will have 3 lagrange coefficients.
[tex]Lg_{1} = \frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})}[/tex]
[tex]Lg_{2} = \frac{(x-x_{1})(x-x_{3})}{(x_{2}-x_{1})(x_{2}-x_{3})}[/tex]
[tex]Lg_{3} = \frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})}[/tex]

Finally your f(x) is given by:
[tex] f(x) = Lg_{1}*y_{1}+Lg_{2}*y_{2}+Lg_{3}*y_{3}[/tex]

Basically the idea is straightforward. For each point you have a corresponding coefficient. So you should loop through all points. For each point's coefficient you need to consider all the other points, so you need to loop through those as well. Basically two for loops will do it. There are ways to optimize your algorithm so that it doesn't perform repeated subtractions.
The following algorithm should give you the general idea (performs in O(n^2) steps):
Code: 
int n //number of points
double xcoords = {x0, x1, x2... x(n-1)}; //the x coordinates of our n points
double ycoords = {y0, y1, y2...y(n-1)}; //the y coordinates of our n points
double xin; //the x whose f(x) we wish to compute
double fx = 0;//the value of f(x)
for(i=0; i<n; i++){
   double Lg = 1;
   for(j=0; j<n; j++){
      if(i != j){
         Lg *= (xin-xcoords[j])/(xcoords[i]-xcoords[j]);
      }
   }
   fx += Lg*ycoords[i];
}
 
Last edited:
-Job- said:
Are you given the points or are you allowed to obtain the points with some strategy? If you have n points, then computing the Lagrange coefficients and interpolating polynomial, f(x), is not very hard. Mostly you should just follow the mathematical expression for the Lagrange coefficients. From what i remember if you have 3 points (x1,y1), (x2,y2), (x3,y3), then you will have 3 lagrange coefficients.
[tex]Lg_{1} = \frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})}[/tex]
[tex]Lg_{2} = \frac{(x-x_{1})(x-x_{3})}{(x_{2}-x_{1})(x_{2}-x_{3})}[/tex]
[tex]Lg_{3} = \frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})}[/tex]

Finally your f(x) is given by:
[tex] f(x) = Lg_{1}*y_{1}+Lg_{2}*y_{2}+Lg_{3}*y_{3}[/tex]

Basically the idea is straightforward. For each point you have a corresponding coefficient. So you should loop through all points. For each point's coefficient you need to consider all the other points, so you need to loop through those as well. Basically two for loops will do it. There are ways to optimize your algorithm so that it doesn't perform repeated subtractions.
The following algorithm should give you the general idea (performs in O(n^2) steps):
Code: 
int n //number of points
double xcoords = {x0, x1, x2... x(n-1)}; //the x coordinates of our n points
double ycoords = {y0, y1, y2...y(n-1)}; //the y coordinates of our n points
double xin; //the x whose f(x) we wish to compute
double fx = 0;//the value of f(x)
for(i=0; i<n; i++){
   double Lg = 1;
   for(j=0; j<n; j++){
      if(i != j){
         Lg *= (xin-xcoords[j])/(xcoords[i]-xcoords[j]);
      }
   }
   fx += Lg*ycoords[i];
}

Thanks for your help. Just to confirm, the program says

Write a program that implemets n-point Lagrange interpolation. Trean n as an input parameter. Apply the program to study the quality of the Lagrange interpolation to function f(x)=sin(x2) initially calculated in 10 uniform points in the interval [0.0, 5.0].

So, does this mean that the user can select the n?
 
Last edited:
Yes, that seems to be the case. You're given the function that you want to model with the interpolation polynomial, [tex]f(x)=sin(x^{2})[/tex]. Apparently you should start with n = 10, and allow the user to enter a larger value for n. Keep in my mind that the value of n includes the endpoints. So for n=10 you have the points x=0, x=5/(10-1), x=2*5/(10-1), x=3*5/(10-1), x=4*5/(10-1), x=5*5/(10-1)... x=(10-1)*5/(10-1)
Basically, the x coordinate of the ith point ([tex]0 \leq i \leq n-1[/tex]) is given by:
[tex]i*\frac{5}{n-1}[/tex]
Then with the x coordinate you get the y coordinate using [tex]f(x)=sin(x^{2})[/tex] and you have your point. I would create a function getIthPoint(i, n) to make your life easier. This way you don't need arrays.
Make sure to go over this in case i made a mistake.
 
Last edited:
-Job- said:
Yes, that seems to be the case. You're given the function that you want to model with the interpolation polynomial, [tex]f(x)=sin(x^{2})[/tex]. Apparently you should start with n = 10, and allow the user to enter a larger value for n. Keep in my mind that the value of n includes the endpoints. So for n=10 you have the points x=0, x=5/(10-1), x=2*5/(10-1), x=3*5/(10-1), x=4*5/(10-1), x=5*5/(10-1)... x=(10-1)*5/(10-1)
Basically, the x coordinate of the ith point ([tex]0 \leq i \leq n-1[/tex]) is given by:
[tex]i*\frac{5}{n-1}[/tex]
Then with the x coordinate you get the y coordinate using [tex]f(x)=sin(x^{2})[/tex] and you have your point. I would create a function getIthPoint(i, n) to make your life easier. This way you don't need arrays.
Make sure to go over this in case i made a mistake.

Hi, I have sort of worked on the code, but I am not getting it work. What am I doing wrong?

//* Lagrange's interpolation */

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>

int main()
{
int n; //Degree of interpolation = (#of nodes -1)
int number;
double* x, *y; //interpolation data
double xx;
double yy; // value to be found = P(xx)
double t; //term in Lagrange's formula

int i,j; //loop indices

count << "Lagrange interpolation program" <<endl;

count << "Enter n input parameter" <<endl;
cin >> n;


if (n<0)
{
count <<"n must be positive!" <<endl;
return(-1);
}
if (n>20)
{
count <<" this program doesn't accept this value" <endl;
return(-1);
}

// count <<"chose the function you want to apply the lagrange interpolation"<<endl;
// cin >> number;


for (i=0; i<=n; i++)
{
x=i*(5/n-1);
y=sin(x*2);


//Calculation

yy=0; //initialization of accumulator
for (i=0; i<=n; i++)
{
xx;

t=y; //initialization of the new term
for (j=0; j<=n; j++)
{
if (j!=i)
{
t*=((xx-x[j])/(x-x[j])); //accumulating product
}
}
yy+=t; //adding the product
}

//Output
count << "result" << xx << yy << endl;

return(0);
}
 
Last edited:
Where do you assign a value to xx(i)? You may want to evaluate you polynomial at many points in your interval (many ~ 100) at the same time you can compute the exact value. Your final print out would be, the interpolated value, the exact value and the difference. This would satisfy the "investigate" part of the problem.

You could also try different interpolation points, both different number and different spacings. (more in fast changing areas, fewer in the slower changing areas etc)
 
From a quick look at your code i see that your f(x) seems to be different. In your code you used:
Code: 
y[i]=sin(x[i]*2);
However, in your initial post you mentioned that your f(x) was:
[tex]f(x)=sin(x^{2})[/tex]
 
-Job- said:
From a quick look at your code i see that your f(x) seems to be different. In your code you used:
Code: 
y[i]=sin(x[i]*2);
However, in your initial post you mentioned that your f(x) was:
[tex]f(x)=sin(x^{2})[/tex]
Sorry, that was a mistake. It is supposed to be sin(x^2), but does everything else seem okay? I am not sure if xx should be 1,2,3... etc...
 
Your xx(i) should be much closer together then your interpolation points. You want to see how well the polynomial interpolation matches the original function, you must compute it at points away from your interpolating points. Take steps in the range of [itex]2^ {-5}[/itex] to [itex]2^{-8}[/itex]
 
Last edited:
  • #10
I did a quick c++ program with the code i posted previously, for testing:

Code: 
#include "stdafx.h"
#include "math.h"
#using <mscorlib.dll>

using namespace System;
double getY(double x){
	return sin(pow(x, 2));
}
double interpolate(double xin, double n){
	double i, j; 
	double fx = 0;
	for(i=0; i<n; i++){
		double Lg = 1;
		double ix = (double)i*5/(n-1);
		double iy = getY(ix);
		for(j=0; j<n; j++){
			if(i!=j){
				double jx = (double)j*5/(n-1);
				Lg *= (xin-jx)/(ix-jx);
			}
		}
		fx += Lg*iy;
	}
	return fx;
}
int _tmain()
{
	double x;
	for(x=0; x<=5; x+=0.05){
		Console::Write("x=");
		Console::Write(x);
		Console::Write("   f(x)=");
		Console::Write(getY(x));
		Console::Write("   P(x)=");
		Console::Write(interpolate(x, 100));
		Console::Write("\n");
	}
	return 0;
}

It's interesting that the above program correctly approximates the value of f(x) with P(x). However there are some values of x where P(x) is way off. I printed the approximations for various x's in the interval [0, 5] and from this data it seems that this behavior is caused by Runge's Phenomenon
Here's the output (notice that P(x) is close to f(x)) for most values in [0, 5]. However in two areas near the endpoints it is way off, and increasing n causes it to be even farther off. This output was printed for n=100)

Code: 
x=0   f(x)=0   P(x)=0
x=0.05   f(x)=0.00249999739583415   P(x)=4.7405971714772
x=0.1   f(x)=0.00999983333416667   P(x)=7.32717058245029
x=0.15   f(x)=0.0224981016105536   P(x)=8.35535776490926
x=0.2   f(x)=0.0399893341866342   P(x)=8.30389975623207
x=0.25   f(x)=0.0624593178423802   P(x)=7.55149792893023
x=0.3   f(x)=0.089878549198011   P(x)=6.39189815050929
x=0.35   f(x)=0.122193852192663   P(x)=5.04733241070067
x=0.4   f(x)=0.159318206614246   P(x)=3.68044203885949
x=0.45   f(x)=0.201118873846073   P(x)=2.40480075249365
x=0.5   f(x)=0.247403959254523   P(x)=1.29415002869104
x=0.55   f(x)=0.297907621896134   P(x)=0.390453673645189
x=0.6   f(x)=0.35227423327509   P(x)=-0.289127018455357
x=0.65   f(x)=0.410041898781764   P(x)=-0.746241477212139
x=0.7   f(x)=0.470625888171158   P(x)=-0.996250617186627
x=0.75   f(x)=0.53330267353602   P(x)=-1.06393441320525
x=0.8   f(x)=0.597195441362392   P(x)=-0.979925507958984
x=0.85   f(x)=0.661262123760472   P(x)=-0.777792890834096
x=0.9   f(x)=0.724287174370143   P(x)=-0.491704375056283
x=0.95   f(x)=0.784878485034067   P(x)=-0.15460115574461
x=1   f(x)=0.841470984807897   P(x)=0.203177845645827
x=1.05   f(x)=0.89233856416221   P(x)=0.554879854762585
x=1.1   f(x)=0.935616001553386   P(x)=0.877984182326532
x=1.15   f(x)=0.969332510867995   P(x)=1.15456004829728
x=1.2   f(x)=0.991458348191687   P(x)=1.37138222153984
x=1.25   f(x)=0.999965585678249   P(x)=1.51984664457424
x=1.3   f(x)=0.992903651094118   P(x)=1.59572458589665
x=1.35   f(x)=0.968489520283355   P(x)=1.59879036789899
x=1.4   f(x)=0.925211520788168   P(x)=1.53235435658385
x=1.45   f(x)=0.86194455517421   P(x)=1.40272967007313
x=1.5   f(x)=0.77807319688792   P(x)=1.21865796634192
x=1.55   f(x)=0.673617587361252   P(x)=0.990716706673673
x=1.6   f(x)=0.549355436427125   P(x)=0.730727460028961
x=1.65   f(x)=0.406931797351062   P(x)=0.451182114847959
x=1.7   f(x)=0.24894678667315   P(x)=0.164701298765956
x=1.75   f(x)=0.0790102167473866   P(x)=-0.116463126687781
x=1.8   f(x)=-0.0982485937451118   P(x)=-0.380871930115067
x=1.85   f(x)=-0.27722754488774   P(x)=-0.618281334851053
x=1.9   f(x)=-0.451465752161427   P(x)=-0.819959441511353
x=1.95   f(x)=-0.613833396107784   P(x)=-0.978932788734592
x=2   f(x)=-0.756802495307931   P(x)=-1.09016136262222
x=2.05   f(x)=-0.872798699517351   P(x)=-1.15064200299717
x=2.1   f(x)=-0.954627771660217   P(x)=-1.15944165959111
x=2.15   f(x)=-0.995962705156195   P(x)=-1.11766332362669
x=2.2   f(x)=-0.991868757310913   P(x)=-1.02834869998609
x=2.25   f(x)=-0.939334638675732   P(x)=-0.896322792250936
x=2.3   f(x)=-0.837769480165098   P(x)=-0.727986547360243
x=2.35   f(x)=-0.689418018619284   P(x)=-0.531064548463635
x=2.4   f(x)=-0.499641883116905   P(x)=-0.31431545374582
x=2.45   f(x)=-0.277014201809769   P(x)=-0.0872134555655052
x=2.5   f(x)=-0.0331792165475613   P(x)=0.140390521812358
x=2.55   f(x)=0.217560782359321   P(x)=0.35861885630802
x=2.6   f(x)=0.458951486377685   P(x)=0.557919285686531
x=2.65   f(x)=0.673781675524276   P(x)=0.729409734488984
x=2.7   f(x)=0.845133411657213   P(x)=0.865206742529471
x=2.75   f(x)=0.957819147348828   P(x)=0.958728546473952
x=2.8   f(x)=0.999902258547975   P(x)=1.00496421787097
x=2.85   f(x)=0.964165036726313   P(x)=1.00070074548939
x=2.9   f(x)=0.849363378505475   P(x)=0.944700566935161
x=2.95   f(x)=0.661095543640036   P(x)=0.837822804267547
x=3   f(x)=0.412118485241771   P(x)=0.683082340715179
x=3.05   f(x)=0.121973473857497   P(x)=0.485641890603811
x=3.1   f(x)=-0.184164779400656   P(x)=0.252733362250751
x=3.15   f(x)=-0.477425197790973   P(x)=-0.00649390614439077
x=3.2   f(x)=-0.727877870349722   P(x)=-0.281199044676839
x=3.25   f(x)=-0.907679875544638   P(x)=-0.559149145508065
x=3.3   f(x)=-0.994432209303193   P(x)=-0.827090135588535
x=3.35   f(x)=-0.97436266119118   P(x)=-1.07118500088798
x=3.4   f(x)=-0.844895943776041   P(x)=-1.27751339125691
x=3.45   f(x)=-0.616170006866291   P(x)=-1.43262455897897
x=3.5   f(x)=-0.311119354981156   P(x)=-1.52413337537844
x=3.55   f(x)=0.0361215260101955   P(x)=-1.54134682851958
x=3.6   f(x)=0.383542755412577   P(x)=-1.47590593107628
x=3.65   f(x)=0.686110730859183   P(x)=-1.32242536085987
x=3.69999999999999   f(x)=0.901675770066375   P(x)=-1.07911041727123
x=3.74999999999999   f(x)=0.997213718805396   P(x)=-0.748328005089993
x=3.79999999999999   f(x)=0.954495430240934   P(x)=-0.337105352528581
x=3.84999999999999   f(x)=0.774208286800343   P(x)=0.142472966637726
x=3.89999999999999   f(x)=0.477637144914053   P(x)=0.673002406801439
x=3.94999999999999   f(x)=0.105267874095546   P(x)=1.23166632243078
x=3.99999999999999   f(x)=-0.287903316665018   P(x)=1.79035397548258
x=4.04999999999999   f(x)=-0.640029556161847   P(x)=2.31603198329057
x=4.09999999999999   f(x)=-0.892129364694352   P(x)=2.77141631217395
x=4.14999999999999   f(x)=-0.998417846377416   P(x)=3.11599585491698
x=4.19999999999999   f(x)=-0.935459140991685   P(x)=3.30746269580768
x=4.24999999999999   f(x)=-0.708278021051184   P(x)=3.30360836509292
x=4.29999999999999   f(x)=-0.35185858693223   P(x)=3.06474971550771
x=4.34999999999999   f(x)=0.0728794083906834   P(x)=2.55675251696862
x=4.39999999999999   f(x)=0.488564765772459   P(x)=1.75472546158502
x=4.44999999999999   f(x)=0.815124497793637   P(x)=0.64746199983705
x=4.49999999999999   f(x)=0.985525111565108   P(x)=-0.757287709903798
x=4.54999999999999   f(x)=0.960459678867108   P(x)=-2.42662746237741
x=4.59999999999999   f(x)=0.738706029703826   P(x)=-4.29531587423681
x=4.64999999999999   f(x)=0.360355131767404   P(x)=-6.25745979870335
x=4.69999999999999   f(x)=-0.0986905140096249   P(x)=-8.1569972322765
x=4.74999999999999   f(x)=-0.540769321526501   P(x)=-9.77686113521475
x=4.79999999999999   f(x)=-0.866851155826269   P(x)=-10.8267099303881
x=4.84999999999999   f(x)=-0.999222150718628   P(x)=-10.9291046553546
x=4.89999999999999   f(x)=-0.901291364088664   P(x)=-9.60400682013398
x=4.94999999999999   f(x)=-0.589339660239126   P(x)=-6.25146496814124
x=4.99999999999999   f(x)=-0.132351750097868   P(x)=-0.132351750099295
Press any key to continue
(notice the round off error in the increment of x, due to 0.05 being a repeating binary sequence, check this thread to see why this is dangerous :smile:)
 
Last edited:
  • #11
-Job- said:
I did a quick c++ program with the code i posted previously, for testing:

Code: 
#include "stdafx.h"
#include "math.h"
#using <mscorlib.dll>

using namespace System;
double getY(double x){
	return sin(pow(x, 2));
}
double interpolate(double xin, double n){
	double i, j; 
	double fx = 0;
	for(i=0; i<n; i++){
		double Lg = 1;
		double ix = (double)i*5/(n-1);
		double iy = getY(ix);
		for(j=0; j<n; j++){
			if(i!=j){
				double jx = (double)j*5/(n-1);
				Lg *= (xin-jx)/(ix-jx);
			}
		}
		fx += Lg*iy;
	}
	return fx;
}
int _tmain()
{
	double x;
	for(x=0; x<=5; x+=0.05){
		Console::Write("x=");
		Console::Write(x);
		Console::Write("   f(x)=");
		Console::Write(getY(x));
		Console::Write("   P(x)=");
		Console::Write(interpolate(x, 100));
		Console::Write("\n");
	}
	return 0;
}

It's interesting that the above program correctly approximates the value of f(x) with P(x). However there are some values of x where P(x) is way off. I printed the approximations for various x's in the interval [0, 5] and from this data it seems that this behavior is caused by Runge's Phenomenon
Here's the output (notice that P(x) is close to f(x)) for most values in [0, 5]. However in two areas near the endpoints it is way off, and increasing n causes it to be even farther off. This output was printed for n=100)

Code: 
x=0   f(x)=0   P(x)=0
x=0.05   f(x)=0.00249999739583415   P(x)=4.7405971714772
x=0.1   f(x)=0.00999983333416667   P(x)=7.32717058245029
x=0.15   f(x)=0.0224981016105536   P(x)=8.35535776490926
x=0.2   f(x)=0.0399893341866342   P(x)=8.30389975623207
x=0.25   f(x)=0.0624593178423802   P(x)=7.55149792893023
x=0.3   f(x)=0.089878549198011   P(x)=6.39189815050929
x=0.35   f(x)=0.122193852192663   P(x)=5.04733241070067
x=0.4   f(x)=0.159318206614246   P(x)=3.68044203885949
x=0.45   f(x)=0.201118873846073   P(x)=2.40480075249365
x=0.5   f(x)=0.247403959254523   P(x)=1.29415002869104
x=0.55   f(x)=0.297907621896134   P(x)=0.390453673645189
x=0.6   f(x)=0.35227423327509   P(x)=-0.289127018455357
x=0.65   f(x)=0.410041898781764   P(x)=-0.746241477212139
x=0.7   f(x)=0.470625888171158   P(x)=-0.996250617186627
x=0.75   f(x)=0.53330267353602   P(x)=-1.06393441320525
x=0.8   f(x)=0.597195441362392   P(x)=-0.979925507958984
x=0.85   f(x)=0.661262123760472   P(x)=-0.777792890834096
x=0.9   f(x)=0.724287174370143   P(x)=-0.491704375056283
x=0.95   f(x)=0.784878485034067   P(x)=-0.15460115574461
x=1   f(x)=0.841470984807897   P(x)=0.203177845645827
x=1.05   f(x)=0.89233856416221   P(x)=0.554879854762585
x=1.1   f(x)=0.935616001553386   P(x)=0.877984182326532
x=1.15   f(x)=0.969332510867995   P(x)=1.15456004829728
x=1.2   f(x)=0.991458348191687   P(x)=1.37138222153984
x=1.25   f(x)=0.999965585678249   P(x)=1.51984664457424
x=1.3   f(x)=0.992903651094118   P(x)=1.59572458589665
x=1.35   f(x)=0.968489520283355   P(x)=1.59879036789899
x=1.4   f(x)=0.925211520788168   P(x)=1.53235435658385
x=1.45   f(x)=0.86194455517421   P(x)=1.40272967007313
x=1.5   f(x)=0.77807319688792   P(x)=1.21865796634192
x=1.55   f(x)=0.673617587361252   P(x)=0.990716706673673
x=1.6   f(x)=0.549355436427125   P(x)=0.730727460028961
x=1.65   f(x)=0.406931797351062   P(x)=0.451182114847959
x=1.7   f(x)=0.24894678667315   P(x)=0.164701298765956
x=1.75   f(x)=0.0790102167473866   P(x)=-0.116463126687781
x=1.8   f(x)=-0.0982485937451118   P(x)=-0.380871930115067
x=1.85   f(x)=-0.27722754488774   P(x)=-0.618281334851053
x=1.9   f(x)=-0.451465752161427   P(x)=-0.819959441511353
x=1.95   f(x)=-0.613833396107784   P(x)=-0.978932788734592
x=2   f(x)=-0.756802495307931   P(x)=-1.09016136262222
x=2.05   f(x)=-0.872798699517351   P(x)=-1.15064200299717
x=2.1   f(x)=-0.954627771660217   P(x)=-1.15944165959111
x=2.15   f(x)=-0.995962705156195   P(x)=-1.11766332362669
x=2.2   f(x)=-0.991868757310913   P(x)=-1.02834869998609
x=2.25   f(x)=-0.939334638675732   P(x)=-0.896322792250936
x=2.3   f(x)=-0.837769480165098   P(x)=-0.727986547360243
x=2.35   f(x)=-0.689418018619284   P(x)=-0.531064548463635
x=2.4   f(x)=-0.499641883116905   P(x)=-0.31431545374582
x=2.45   f(x)=-0.277014201809769   P(x)=-0.0872134555655052
x=2.5   f(x)=-0.0331792165475613   P(x)=0.140390521812358
x=2.55   f(x)=0.217560782359321   P(x)=0.35861885630802
x=2.6   f(x)=0.458951486377685   P(x)=0.557919285686531
x=2.65   f(x)=0.673781675524276   P(x)=0.729409734488984
x=2.7   f(x)=0.845133411657213   P(x)=0.865206742529471
x=2.75   f(x)=0.957819147348828   P(x)=0.958728546473952
x=2.8   f(x)=0.999902258547975   P(x)=1.00496421787097
x=2.85   f(x)=0.964165036726313   P(x)=1.00070074548939
x=2.9   f(x)=0.849363378505475   P(x)=0.944700566935161
x=2.95   f(x)=0.661095543640036   P(x)=0.837822804267547
x=3   f(x)=0.412118485241771   P(x)=0.683082340715179
x=3.05   f(x)=0.121973473857497   P(x)=0.485641890603811
x=3.1   f(x)=-0.184164779400656   P(x)=0.252733362250751
x=3.15   f(x)=-0.477425197790973   P(x)=-0.00649390614439077
x=3.2   f(x)=-0.727877870349722   P(x)=-0.281199044676839
x=3.25   f(x)=-0.907679875544638   P(x)=-0.559149145508065
x=3.3   f(x)=-0.994432209303193   P(x)=-0.827090135588535
x=3.35   f(x)=-0.97436266119118   P(x)=-1.07118500088798
x=3.4   f(x)=-0.844895943776041   P(x)=-1.27751339125691
x=3.45   f(x)=-0.616170006866291   P(x)=-1.43262455897897
x=3.5   f(x)=-0.311119354981156   P(x)=-1.52413337537844
x=3.55   f(x)=0.0361215260101955   P(x)=-1.54134682851958
x=3.6   f(x)=0.383542755412577   P(x)=-1.47590593107628
x=3.65   f(x)=0.686110730859183   P(x)=-1.32242536085987
x=3.69999999999999   f(x)=0.901675770066375   P(x)=-1.07911041727123
x=3.74999999999999   f(x)=0.997213718805396   P(x)=-0.748328005089993
x=3.79999999999999   f(x)=0.954495430240934   P(x)=-0.337105352528581
x=3.84999999999999   f(x)=0.774208286800343   P(x)=0.142472966637726
x=3.89999999999999   f(x)=0.477637144914053   P(x)=0.673002406801439
x=3.94999999999999   f(x)=0.105267874095546   P(x)=1.23166632243078
x=3.99999999999999   f(x)=-0.287903316665018   P(x)=1.79035397548258
x=4.04999999999999   f(x)=-0.640029556161847   P(x)=2.31603198329057
x=4.09999999999999   f(x)=-0.892129364694352   P(x)=2.77141631217395
x=4.14999999999999   f(x)=-0.998417846377416   P(x)=3.11599585491698
x=4.19999999999999   f(x)=-0.935459140991685   P(x)=3.30746269580768
x=4.24999999999999   f(x)=-0.708278021051184   P(x)=3.30360836509292
x=4.29999999999999   f(x)=-0.35185858693223   P(x)=3.06474971550771
x=4.34999999999999   f(x)=0.0728794083906834   P(x)=2.55675251696862
x=4.39999999999999   f(x)=0.488564765772459   P(x)=1.75472546158502
x=4.44999999999999   f(x)=0.815124497793637   P(x)=0.64746199983705
x=4.49999999999999   f(x)=0.985525111565108   P(x)=-0.757287709903798
x=4.54999999999999   f(x)=0.960459678867108   P(x)=-2.42662746237741
x=4.59999999999999   f(x)=0.738706029703826   P(x)=-4.29531587423681
x=4.64999999999999   f(x)=0.360355131767404   P(x)=-6.25745979870335
x=4.69999999999999   f(x)=-0.0986905140096249   P(x)=-8.1569972322765
x=4.74999999999999   f(x)=-0.540769321526501   P(x)=-9.77686113521475
x=4.79999999999999   f(x)=-0.866851155826269   P(x)=-10.8267099303881
x=4.84999999999999   f(x)=-0.999222150718628   P(x)=-10.9291046553546
x=4.89999999999999   f(x)=-0.901291364088664   P(x)=-9.60400682013398
x=4.94999999999999   f(x)=-0.589339660239126   P(x)=-6.25146496814124
x=4.99999999999999   f(x)=-0.132351750097868   P(x)=-0.132351750099295
Press any key to continue
(notice the round off error in the increment of x, due to 0.05 being a repeating binary sequence, check this thread to see why this is dangerous :smile:)

Thanks, that really cleared up a lot of stuff. Yeah, I have used arrays in my code, and it was getting too complicated, but I am still not sure why my code didn't work with the arrays.

Now, if I were to compare lagrange interpolation, to say, cubic spline interpolation, then it's be better to make a new function, right? Also, the x and y values are going to be pretty much the ones from the sinx^2 function, right? I am hoping to somehow use the lagrange interpolation function to sort of minimize the code I have to write for the cubic spline.
 
Last edited:
  • #12
I just want to correct something. The output that i posted was actually for n=10, that's why P(x) and f(x) are somewhat off. I meant to post the output with n=100 which is much better and you can clearly see Runge's Phenomenon:

Code: 
x=0   f(x)=0   P(x)=0
x=0.05   f(x)=0.00249999739583415   P(x)=-60907171.41443
x=0.1   f(x)=0.00999983333416667   P(x)=-38762.241528455
x=0.15   f(x)=0.0224981016105536   P(x)=133888.976920124
x=0.2   f(x)=0.0399893341866342   P(x)=-47477.1557530205
x=0.25   f(x)=0.0624593178423802   P(x)=-676.321658577803
x=0.3   f(x)=0.089878549198011   P(x)=-70.0563415465577
x=0.35   f(x)=0.122193852192663   P(x)=-10.3323357821072
x=0.4   f(x)=0.159318206614246   P(x)=-0.743757999126726
x=0.45   f(x)=0.201118873846073   P(x)=0.242422729546865
x=0.5   f(x)=0.247403959254523   P(x)=0.234926987631692
x=0.55   f(x)=0.297907621896134   P(x)=0.297937142509013
x=0.6   f(x)=0.35227423327509   P(x)=0.352142785624735
x=0.65   f(x)=0.410041898781764   P(x)=0.410028523227198
x=0.7   f(x)=0.470625888171158   P(x)=0.470635245914266
x=0.75   f(x)=0.53330267353602   P(x)=0.533302801399181
x=0.8   f(x)=0.597195441362392   P(x)=0.59719531855193
x=0.85   f(x)=0.661262123760472   P(x)=0.661262069516314
x=0.9   f(x)=0.724287174370143   P(x)=0.724287171972466
x=0.95   f(x)=0.784878485034067   P(x)=0.784878488139338
x=1   f(x)=0.841470984807897   P(x)=0.84147098533634
x=1.05   f(x)=0.89233856416221   P(x)=0.892338563879778
x=1.1   f(x)=0.935616001553386   P(x)=0.935616001445517
x=1.15   f(x)=0.969332510867995   P(x)=0.969332510854974
x=1.2   f(x)=0.991458348191687   P(x)=0.99145834818017
x=1.25   f(x)=0.999965585678249   P(x)=0.999965585674027
x=1.3   f(x)=0.992903651094118   P(x)=0.992903651094568
x=1.35   f(x)=0.968489520283355   P(x)=0.968489520283317
x=1.4   f(x)=0.925211520788168   P(x)=0.925211520787792
x=1.45   f(x)=0.86194455517421   P(x)=0.861944555174197
x=1.5   f(x)=0.77807319688792   P(x)=0.778073196887921
x=1.55   f(x)=0.673617587361252   P(x)=0.67361758736126
x=1.6   f(x)=0.549355436427125   P(x)=0.549355436427112
x=1.65   f(x)=0.406931797351062   P(x)=0.406931797351064
x=1.7   f(x)=0.24894678667315   P(x)=0.248946786673147
x=1.75   f(x)=0.0790102167473866   P(x)=0.0790102167473904
x=1.8   f(x)=-0.0982485937451118   P(x)=-0.0982485937451123
x=1.85   f(x)=-0.27722754488774   P(x)=-0.27722754488774
x=1.9   f(x)=-0.451465752161427   P(x)=-0.451465752161426
x=1.95   f(x)=-0.613833396107784   P(x)=-0.613833396107783
x=2   f(x)=-0.756802495307931   P(x)=-0.75680249530793
x=2.05   f(x)=-0.872798699517351   P(x)=-0.872798699517349
x=2.1   f(x)=-0.954627771660217   P(x)=-0.954627771660218
x=2.15   f(x)=-0.995962705156195   P(x)=-0.995962705156194
x=2.2   f(x)=-0.991868757310913   P(x)=-0.991868757310912
x=2.25   f(x)=-0.939334638675732   P(x)=-0.939334638675732
x=2.3   f(x)=-0.837769480165098   P(x)=-0.837769480165098
x=2.35   f(x)=-0.689418018619284   P(x)=-0.689418018619284
x=2.4   f(x)=-0.499641883116905   P(x)=-0.499641883116904
x=2.45   f(x)=-0.277014201809769   P(x)=-0.277014201809769
x=2.5   f(x)=-0.0331792165475613   P(x)=-0.0331792165475615
x=2.55   f(x)=0.217560782359321   P(x)=0.217560782359322
x=2.6   f(x)=0.458951486377685   P(x)=0.458951486377685
x=2.65   f(x)=0.673781675524276   P(x)=0.673781675524275
x=2.7   f(x)=0.845133411657213   P(x)=0.845133411657214
x=2.75   f(x)=0.957819147348828   P(x)=0.957819147348828
x=2.8   f(x)=0.999902258547975   P(x)=0.999902258547975
x=2.85   f(x)=0.964165036726313   P(x)=0.964165036726314
x=2.9   f(x)=0.849363378505475   P(x)=0.849363378505476
x=2.95   f(x)=0.661095543640036   P(x)=0.661095543640036
x=3   f(x)=0.412118485241771   P(x)=0.412118485241773
x=3.05   f(x)=0.121973473857497   P(x)=0.121973473857497
x=3.1   f(x)=-0.184164779400656   P(x)=-0.184164779400655
x=3.15   f(x)=-0.477425197790973   P(x)=-0.477425197790976
x=3.2   f(x)=-0.727877870349722   P(x)=-0.727877870349721
x=3.25   f(x)=-0.907679875544638   P(x)=-0.907679875544638
x=3.3   f(x)=-0.994432209303193   P(x)=-0.994432209303187
x=3.35   f(x)=-0.97436266119118   P(x)=-0.97436266119119
x=3.4   f(x)=-0.844895943776041   P(x)=-0.844895943776045
x=3.45   f(x)=-0.616170006866291   P(x)=-0.616170006866324
x=3.5   f(x)=-0.311119354981156   P(x)=-0.311119354981138
x=3.55   f(x)=0.0361215260101955   P(x)=0.0361215260101664
x=3.6   f(x)=0.383542755412577   P(x)=0.383542755412759
x=3.65   f(x)=0.686110730859183   P(x)=0.686110730857785
x=3.69999999999999   f(x)=0.901675770066375   P(x)=0.901675770068809
x=3.74999999999999   f(x)=0.997213718805396   P(x)=0.997213718809361
x=3.79999999999999   f(x)=0.954495430240934   P(x)=0.954495430259517
x=3.84999999999999   f(x)=0.774208286800343   P(x)=0.77420828683213
x=3.89999999999999   f(x)=0.477637144914053   P(x)=0.477637144984646
x=3.94999999999999   f(x)=0.105267874095546   P(x)=0.105267874149488
x=3.99999999999999   f(x)=-0.287903316665018   P(x)=-0.287903314886773
x=4.04999999999999   f(x)=-0.640029556161847   P(x)=-0.64002955654917
x=4.09999999999999   f(x)=-0.892129364694352   P(x)=-0.892129380112669
x=4.14999999999999   f(x)=-0.998417846377416   P(x)=-0.998417955099416
x=4.19999999999999   f(x)=-0.935459140991685   P(x)=-0.935459246420381
x=4.24999999999999   f(x)=-0.708278021051184   P(x)=-0.708279101677588
x=4.29999999999999   f(x)=-0.35185858693223   P(x)=-0.351858017815286
x=4.34999999999999   f(x)=0.0728794083906834   P(x)=0.0728393961814019
x=4.39999999999999   f(x)=0.488564765772459   P(x)=0.488707570655531
x=4.44999999999999   f(x)=0.815124497793637   P(x)=0.814437266123537
x=4.49999999999999   f(x)=0.985525111565108   P(x)=1.00134556766263
x=4.54999999999999   f(x)=0.960459678867108   P(x)=1.00535734265019
x=4.59999999999999   f(x)=0.738706029703826   P(x)=-0.916344799364983
x=4.64999999999999   f(x)=0.360355131767404   P(x)=5.55297239198827
x=4.69999999999999   f(x)=-0.0986905140096249   P(x)=-19.3321227061742
x=4.74999999999999   f(x)=-0.540769321526501   P(x)=-1764.86886307135
x=4.79999999999999   f(x)=-0.866851155826269   P(x)=-2493.22965470624
x=4.84999999999999   f(x)=-0.999222150718628   P(x)=-218569.956219739
x=4.89999999999999   f(x)=-0.901291364088664   P(x)=8161924.37235928
x=4.94999999999999   f(x)=-0.589339660239126   P(x)=179083586.20583
x=4.99999999999999   f(x)=-0.132351750097868   P(x)=-0.320717299305203
Press any key to continue

For cubic splines i suppose you can use the same strategy for obtaining the points.
 
Last edited:
  • #13
Nice work Job!

Include a difference column (F(x) - PN) This will revel some structure in the interpolation errors.

Yes, you should stay away from the end points of High order polynoimials they can get pretty snaky. Higher order is not always better, a lot depends on the behavior of your function in the region of your approximation. The end effects may cause some trouble linking to a cubic spline. I am not real sure what the OP has in mind so can't say much more without more info.

You want to use the http://www.cse.uiuc.edu/eot/modules/interpolation/chbshvp/" points as your interpolating set. This is the minmial error.

Happy computing.:smile:
 
Last edited by a moderator:
  • #14
Does anyone by any chance have that in fortran? i will be needing it and i unfortunately know almost nothing about fortran whatsoever
 
  • #15
Lagrange Interpolation using Newton

Hello
Can you help me with an inplementation of lagrange interpolating polynom using the Newton's formula? (that with divided differences). or smth appropiate to that, i want to calculate the lagrange interp. polynom, but not using the barycentric formula
Thanks :)

p.s. the Input is the nodes to which i want to calculate the polynom
 

Similar threads

How Can Lagrange Interpolation Be Implemented in Software Development?
  • · Replies 1 ·
Replies
1
Views
7K
Fortran Interpolating & Integrating Data from File to 2D Array
  • · Replies 10 ·
Replies
10
Views
26K
C/C++ How to Implement Polynomial Operations in C++?
  • · Replies 5 ·
Replies
5
Views
3K
Fortran Fortran, interpolating polynomial program doesn't compile
  • · Replies 25 ·
Replies
25
Views
10K
C/C++ Puzzle with C++ constructor and overloaded operator
  • · Replies 23 ·
Replies
23
Views
3K
Fortran Simulating Contact Force in a Moving Dashpot System Using FORTRAN
  • · Replies 26 ·
Replies
26
Views
4K
Solving Newton interpolation with c
  • · Replies 8 ·
Replies
8
Views
8K
Comp Sci [Fortran] making a more efficient bilinear interpolation
  • · Replies 11 ·
Replies
11
Views
6K
MHB Maximum error for the Lagrange interpolating polynomials
  • · Replies 8 ·
Replies
8
Views
2K
Fortran Help with Fortran 77 Programming for Two-Dimensional Axisymmetric Flow Study
  • · Replies 6 ·
Replies
6
Views
4K