Lagrange interpolation over a larger set of data 1

dibas · Sep 25, 2012

Hi all. I have a Lagrange program that uses a 4-point data set to interpolate some non-tabulated value
in between any of those data points. I want to extend this script to act on a large input file, call
it file1.dat, where the data is arranged in two columns Xin Yin. Then the script should interpolate
for a range of set values Xnew to get Ynew (i.e. a new output data file file2.dat with two columns
Xnew Ynew). I need to use this to investigate the convergence of some simulations which I run for two
different resolutions whereby the input gridpoints do not coincide for those two resolutions. So, what
I want to do is to interpolate the higher resolution data (file1.dat) onto the lower resolution one so
that it can be possible to plot some quantity from the two different runs on the same graph and study
clearly its convergence.

Code:

program Lagrange_int

! A program that approximates a value of a function/experimental quantity at a point based
! on a given set of data

  implicit none
  integer, parameter :: n = 4 ! number of data points:length of each data vector
  integer :: i, j
  double precision, dimension(n) :: xvect, yvect! vectors of known data: x,y-values
  double precision, dimension(n) :: L  ! Lagrange cardinal function
  double precision :: z, Pz  ! the set evaluation point "z" and the value of the 
                             ! ploynomial/function/quantity at that point

! Initializations of z, Pz and L
  z = 3.0d0  ! the point of evaluation
  Pz = 0.0d0 ! initializing the polynomia value at z
  L  = 1.0d0  ! initalizing the vector of cardinal functions to 1

! The data
  data (xvect(i), i=1,n)/3.2d0,2.7d0,1.0d0,4.8d0/
  data (yvect(i), i=1,n)/22.d0,17.8d0,14.2d0,38.3d0/

! Performing the interpolation
  do i = 1,n
     do j = 1, n
        if (i /= j) then
           L(i) = ( (z - xvect(j)) / (xvect(i) - xvect(j)) )* L(i) ! part of L(i)
        end if
     end do
     Pz = Pz + L(i)*yvect(i) ! update Pz ~ f(z)
  end do

  write(*,*) "The approximate value of y(",z,")=", Pz

end program Lagrange_int

gummibaer · Sep 25, 2012

And what is your problem related to Fortran?

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

mikrom · Sep 25, 2012

Hi dibas,
if I good understand, for your given x you want to find four x-values
x[i-1], x, x[i+1], x[i+2] with corresponding y-values y[i-1], y, y[i+1], y[i+2], so that:
x[i-1] < x < x < x[i+1] < x[i+2]
and compute y(x) using the 3-rd order interpolation polynomial.
But IMO for this reason your values have to be ordered.
In your example above you have
data (xvect(i), i=1,n)/3.2d0,2.7d0,1.0d0,4.8d0/
data (yvect(i), i=1,n)/22.d0,17.8d0,14.2d0,38.3d0/
but IMO it should be ordered like this
data (xvect(i), i=1,n)/1.0d0,2.7d0,3.2d0,4.8d0/
data (yvect(i), i=1,n)/14.2d0,17.8d0,22.d0,38.3d0/

Please post an example of your data.

mikrom · Sep 26, 2012

Hi dibas,
1. You have to wrap your 4-point computation of interpolation polynomial of 3rd order into a procedure.
2. For a given point z, find it's position idx in x-array, so that x[idx] <= z <= x[idx+1].
3. Take the 4 adjacent points from x-array and y-array, e.g. with indices idx-1, idx, idx+1, idx+2 and compute the interpolated value y(z) using the procedure you created in 1.step

dibas · Sep 26, 2012

Hi mikrom,

I am not able to attach the files themselves. So I'll just paste the data below.
The file.dat.400 on is the higher resolution
which is the main input file (I call the y-column there Y_hi). The file.dat.200
is the lower resolution file from which I only need the X's. The aim is to
have the f90 script interpolate Y_hi on the range of points X_low. The X's in
both files are in ascending order.

file.dat.400
0.0250000000000000 -0.0000000000021316
0.1250000000000000 -0.0000000000073896
0.2250000000000000 -0.0000000000280171
0.3250000000000000 -0.0000000000934364
0.4250000000000000 -0.0000000002827281
0.5250000000000000 -0.0000000007933759
0.6250000000000000 -0.0000000020969573
0.7250000000000001 -0.0000000052747691
0.8250000000000001 -0.0000000127158297
0.9250000000000000 -0.0000000295274628
1.0250000000000001 -0.0000000662901387
1.1250000000000000 -0.0000001442749466
1.2250000000000001 -0.0000003050139968
1.3250000000000002 -0.0000006273103797
1.4250000000000000 -0.0000012565028316
1.5250000000000001 -0.0000024531579726
1.6250000000000000 -0.0000046713149151
1.7250000000000001 -0.0000086796764790
1.8250000000000002 -0.0000157420695097
1.9250000000000000 -0.0000278748093975
2.0249999999999999 -0.0000481963644928
2.1250000000000000 -0.0000813762167883
2.2250000000000001 -0.0001341719687054
2.3250000000000002 -0.0002160136869760
2.4250000000000003 -0.0003395508217744
2.5250000000000004 -0.0005210215467427
2.6250000000000000 -0.0007802438047672
2.7250000000000001 -0.0011399750473707
2.8250000000000002 -0.0016243637915730
2.9250000000000003 -0.0022562458481933
3.0250000000000004 -0.0030531466405447
3.1250000000000000 -0.0040220557957028
3.2250000000000001 -0.0051533404691712
3.3250000000000002 -0.0064145311244965
3.4250000000000003 -0.0077450854788916
3.5250000000000004 -0.0090535193797976
3.6250000000000000 -0.0102183760674242
3.7250000000000001 -0.0110942840803917
3.8250000000000002 -0.0115237672978878
3.9250000000000003 -0.0113545326473761
4.0250000000000004 -0.0104607853093806
4.1250000000000000 -0.0087659203649287
4.2250000000000005 -0.0062629951598457
4.3250000000000002 -0.0030289881375843
4.4249999999999998 0.0007707827486166
4.5250000000000004 0.0048905031726057
4.6250000000000000 0.0090244690514987
4.7250000000000005 0.0128386505656734
4.8250000000000002 0.0160091798191104
4.9250000000000007 0.0182620896621704
5.0250000000000004 0.0194080003009810
5.1250000000000000 0.0193660385626524
5.2250000000000005 0.0181730163146236
5.3250000000000002 0.0159764492249323
5.4250000000000007 0.0130128275078594
5.5250000000000004 0.0095750550000660
5.6250000000000000 0.0059746360026926
5.7250000000000005 0.0025047053287008
5.8250000000000002 -0.0005906645743532
5.9250000000000007 -0.0031370507621734
6.0250000000000004 -0.0050384609534960
6.1250000000000000 -0.0062747718344226
6.2250000000000005 -0.0068900523576977
6.3250000000000002 -0.0069750806425784
6.4250000000000007 -0.0066476658640141
6.5250000000000004 -0.0060339978104405
6.6250000000000000 -0.0052533805599960
6.7250000000000005 -0.0044076308833498
6.8250000000000002 -0.0035753871607535
6.9250000000000007 -0.0028107630154732
7.0250000000000004 -0.0021452824606435
7.1250000000000000 -0.0015918539693506
7.2250000000000005 -0.0011496194508172
7.3250000000000002 -0.0008087579359079
7.4250000000000007 -0.0005546369059065
7.5250000000000004 -0.0003710084635589
7.6250000000000000 -0.0002421928699879
7.7250000000000005 -0.0001543569380895
7.8250000000000002 -0.0000960807943772
7.9250000000000007 -0.0000584291421491
8.0250000000000004 -0.0000347235489245
8.1250000000000000 -0.0000201708802429
8.2249999999999996 -0.0000114557596824
8.3250000000000011 -0.0000063621363542
8.4250000000000007 -0.0000034556914339
8.5250000000000004 -0.0000018360523927
8.6250000000000000 -0.0000009543596715
8.7249999999999996 -0.0000004853655049
8.8250000000000011 -0.0000002415480285
8.9250000000000007 -0.0000001176414302
9.0250000000000004 -0.0000000560762295
9.1250000000000000 -0.0000000261634362
9.2249999999999996 -0.0000000119492964
9.3250000000000011 -0.0000000053426218
9.4250000000000007 -0.0000000023385542
9.5250000000000004 -0.0000000010022513
9.6250000000000000 -0.0000000004205977
9.7250000000000014 -0.0000000001728283
9.8250000000000011 -0.0000000000694922
9.9250000000000007 -0.0000000000274704
10.0250000000000004 -0.0000000000106125
10.1250000000000000 -0.0000000000039148
10.2250000000000014 -0.0000000000014808
10.3250000000000011 -0.0000000000005226
10.4250000000000007 -0.0000000000002626
10.5250000000000004 -0.0000000000000876
10.6250000000000000 0.0000000000000892
10.7250000000000014 0.0000000000000000
10.8250000000000011 0.0000000000000000
10.9250000000000007 0.0000000000000000
11.0250000000000004 0.0000000000000000
11.1250000000000000 0.0000000000000000
11.2250000000000014 0.0000000000000000
11.3250000000000011 0.0000000000000000
11.4250000000000007 0.0000000000000000
11.5250000000000004 0.0000000000000000
11.6250000000000000 0.0000000000000000
11.7250000000000014 0.0000000000000000
11.8250000000000011 0.0000000000000000
11.9250000000000007 0.0000000000000000
12.0250000000000004 0.0000000000000000
12.1250000000000000 0.0000000000000000
12.2250000000000014 0.0000000000000000
12.3250000000000011 0.0000000000000000
12.4250000000000007 0.0000000000000000
12.5250000000000004 0.0000000000000000
12.6250000000000000 0.0000000000000000
12.7250000000000014 0.0000000000000000
12.8250000000000011 0.0000000000000000
12.9250000000000007 0.0000000000000000
13.0250000000000004 0.0000000000000000
13.1250000000000000 0.0000000000000000
13.2250000000000014 0.0000000000000000
13.3250000000000011 0.0000000000000000
13.4250000000000007 0.0000000000000000
13.5250000000000004 0.0000000000000000
13.6250000000000000 0.0000000000000000
13.7250000000000014 0.0000000000000000
13.8250000000000011 0.0000000000000000
13.9250000000000007 0.0000000000000000
14.0250000000000004 0.0000000000000000
14.1250000000000000 0.0000000000000000
14.2250000000000014 0.0000000000000000
14.3250000000000011 0.0000000000000000
14.4250000000000007 0.0000000000000000
14.5250000000000004 0.0000000000000000
14.6250000000000000 0.0000000000000000
14.7250000000000014 0.0000000000000000
14.8250000000000011 0.0000000000000000
14.9250000000000007 0.0000000000000000
15.0250000000000004 0.0000000000000000
15.1250000000000000 0.0000000000000000
15.2250000000000014 0.0000000000000000
15.3250000000000011 0.0000000000000000
15.4250000000000007 0.0000000000000000
15.5250000000000004 0.0000000000000000
15.6250000000000000 0.0000000000000000
15.7250000000000014 0.0000000000000000
15.8250000000000011 0.0000000000000000
15.9250000000000007 0.0000000000000000
16.0250000000000021 0.0000000000000000
16.1250000000000000 0.0000000000000000
16.2250000000000014 0.0000000000000000
16.3249999999999993 0.0000000000000000
16.4250000000000007 0.0000000000000000
16.5250000000000021 0.0000000000000000
16.6250000000000000 0.0000000000000000
16.7250000000000014 0.0000000000000000
16.8249999999999993 0.0000000000000000
16.9250000000000007 0.0000000000000000
17.0250000000000021 0.0000000000000000
17.1250000000000000 0.0000000000000000
17.2250000000000014 0.0000000000000000
17.3249999999999993 0.0000000000000000
17.4250000000000007 0.0000000000000000
17.5250000000000021 0.0000000000000000
17.6250000000000000 0.0000000000000000
17.7250000000000014 0.0000000000000000
17.8249999999999993 0.0000000000000000
17.9250000000000007 0.0000000000000000
18.0250000000000021 0.0000000000000000
18.1250000000000000 0.0000000000000000
18.2250000000000014 0.0000000000000000
18.3249999999999993 0.0000000000000000
18.4250000000000007 0.0000000000000000
18.5250000000000021 0.0000000000000000
18.6250000000000000 0.0000000000000000
18.7250000000000014 0.0000000000000000
18.8249999999999993 0.0000000000000000
18.9250000000000007 0.0000000000000000
19.0250000000000021 0.0000000000000000
19.1250000000000000 0.0000000000000000
19.2250000000000014 0.0000000000000000
19.3250000000000028 0.0000000000000000
19.4250000000000007 0.0000000000000000
19.5250000000000021 0.0000000000000000
19.6250000000000000 0.0000000000000000
19.7250000000000014 0.0000000000000000
19.8250000000000028 0.0000000000000000
19.9250000000000007 0.0000000000000000
20.0250000000000021 -0.0000000000000000

file.dat.200
0.0500000000000000 -0.0000000000039302
0.1500000000000000 -0.0000000000129230
0.2500000000000000 -0.0000000000442446
0.3500000000000000 -0.0000000001388984
0.4500000000000000 -0.0000000004034995
0.5500000000000000 -0.0000000010993953
0.6500000000000000 -0.0000000028400991
0.7500000000000000 -0.0000000070115247
0.8500000000000001 -0.0000000166373569
0.9500000000000001 -0.0000000381050127
1.0500000000000000 -0.0000000845017960
1.1500000000000001 -0.0000001818617624
1.2500000000000000 -0.0000003805048800
1.3500000000000001 -0.0000007749720717
1.4500000000000002 -0.0000015379500907
1.5500000000000000 -0.0000029760914908
1.6500000000000001 -0.0000056186809939
1.7500000000000000 -0.0000103533147776
1.8500000000000001 -0.0000186253247534
1.9500000000000002 -0.0000327181733159
2.0500000000000003 -0.0000561283051311
2.1499999999999999 -0.0000940371788366
2.2500000000000000 -0.0001538623287943
2.3500000000000001 -0.0002458358004588
2.4500000000000002 -0.0003835115180412
2.5500000000000003 -0.0005840460360961
2.6500000000000004 -0.0008680381206984
2.7500000000000000 -0.0012586667671867
2.8500000000000001 -0.0017798561539483
2.9500000000000002 -0.0024532450410596
3.0500000000000003 -0.0032938707923701
3.1500000000000004 -0.0043047076580638
3.2500000000000000 -0.0054705177089052
3.3500000000000001 -0.0067518438504252
3.4500000000000002 -0.0080803271702699
3.5500000000000003 -0.0093567666626049
3.6500000000000004 -0.0104533477128409
3.7500000000000000 -0.0112211525703447
3.8500000000000001 -0.0115033880420398
3.9500000000000002 -0.0111537619629799
4.0499999999999998 -0.0100582497365971
4.1500000000000004 -0.0081573469593807
4.2500000000000000 -0.0054650897022350
4.3500000000000005 -0.0020809151411711
4.4500000000000002 0.0018089836628689
4.5499999999999998 0.0059427607206173
4.6500000000000004 0.0100063987389164
4.7500000000000000 0.0136671867727905
4.8500000000000005 0.0166127180814471
4.9500000000000002 0.0185895394701183
5.0500000000000007 0.0194352599212583
5.1500000000000004 0.0190988142247208
5.2500000000000000 0.0176455397042999
5.3500000000000005 0.0152463671630304
5.4500000000000002 0.0121531926558088
5.5500000000000007 0.0086648004742898
5.6500000000000004 0.0050890844814179
5.7500000000000000 0.0017075282895704
5.8500000000000005 -0.0012530070195867
5.9500000000000002 -0.0036378913706743
6.0500000000000007 -0.0053705857034563
6.1500000000000004 -0.0064476818817114
6.2500000000000000 -0.0069256473933564
6.3500000000000005 -0.0069026826671032
6.4500000000000002 -0.0064992211066456
6.5500000000000007 -0.0058400844792436
6.6500000000000004 -0.0050403789806272
6.7500000000000000 -0.0041961454852472
6.8500000000000005 -0.0033797954756598
6.9500000000000002 -0.0026396311597496
7.0500000000000007 -0.0020023326550469
7.1500000000000004 -0.0014771836066909
7.2500000000000000 -0.0010609299454384
7.3500000000000005 -0.0007424304183129
7.4500000000000002 -0.0005065701962238
7.5500000000000007 -0.0003371991049358
7.6500000000000004 -0.0002190811614687
7.7500000000000000 -0.0001389867207503
7.8500000000000005 -0.0000861279158221
7.9500000000000002 -0.0000521493788262
8.0500000000000007 -0.0000308605903351
8.1500000000000004 -0.0000178529689753
8.2500000000000000 -0.0000100985080672
8.3499999999999996 -0.0000055862958691
8.4500000000000011 -0.0000030226055029
8.5500000000000007 -0.0000015999018299
8.6500000000000004 -0.0000008285462121
8.7500000000000000 -0.0000004198589158
8.8499999999999996 -0.0000002082093006
8.9500000000000011 -0.0000001010530628
9.0500000000000007 -0.0000000480054913
9.1500000000000004 -0.0000000223233525
9.2500000000000000 -0.0000000101622120
9.3499999999999996 -0.0000000045290538
9.4500000000000011 -0.0000000019762331
9.5500000000000007 -0.0000000008443848
9.6500000000000004 -0.0000000003532082
9.7500000000000000 -0.0000000001447481
9.8500000000000014 -0.0000000000580297
9.9500000000000011 -0.0000000000228373
10.0500000000000007 -0.0000000000087533
10.1500000000000004 -0.0000000000033257
10.2500000000000000 -0.0000000000012179
10.3500000000000014 -0.0000000000004351
10.4500000000000011 -0.0000000000001522
10.5500000000000007 -0.0000000000000656
10.6500000000000004 0.0000000000000004
10.7500000000000000 -0.0000000000000220
10.8500000000000014 0.0000000000000000
10.9500000000000011 0.0000000000000000
11.0500000000000007 0.0000000000000000
11.1500000000000004 0.0000000000000000
11.2500000000000000 0.0000000000000000
11.3500000000000014 0.0000000000000000
11.4500000000000011 0.0000000000000000
11.5500000000000007 0.0000000000000000
11.6500000000000004 0.0000000000000000
11.7500000000000000 0.0000000000000000
11.8500000000000014 0.0000000000000000
11.9500000000000011 0.0000000000000000
12.0500000000000007 0.0000000000000000
12.1500000000000004 0.0000000000000000
12.2500000000000000 0.0000000000000000
12.3500000000000014 0.0000000000000000
12.4500000000000011 0.0000000000000000
12.5500000000000007 0.0000000000000000
12.6500000000000004 0.0000000000000000
12.7500000000000000 0.0000000000000000
12.8500000000000014 0.0000000000000000
12.9500000000000011 0.0000000000000000
13.0500000000000007 0.0000000000000000
13.1500000000000004 0.0000000000000000
13.2500000000000000 0.0000000000000000
13.3500000000000014 0.0000000000000000
13.4500000000000011 0.0000000000000000
13.5500000000000007 0.0000000000000000
13.6500000000000004 0.0000000000000000
13.7500000000000000 0.0000000000000000
13.8500000000000014 0.0000000000000000
13.9500000000000011 0.0000000000000000
14.0500000000000007 0.0000000000000000
14.1500000000000004 0.0000000000000000
14.2500000000000000 0.0000000000000000
14.3500000000000014 0.0000000000000000
14.4500000000000011 0.0000000000000000
14.5500000000000007 0.0000000000000000
14.6500000000000004 0.0000000000000000
14.7500000000000000 0.0000000000000000
14.8500000000000014 0.0000000000000000
14.9500000000000011 0.0000000000000000
15.0500000000000007 0.0000000000000000
15.1500000000000004 0.0000000000000000
15.2500000000000000 0.0000000000000000
15.3500000000000014 0.0000000000000000
15.4500000000000011 0.0000000000000000
15.5500000000000007 0.0000000000000000
15.6500000000000004 0.0000000000000000
15.7500000000000000 0.0000000000000000
15.8500000000000014 0.0000000000000000
15.9500000000000011 0.0000000000000000
16.0500000000000007 0.0000000000000000
16.1500000000000021 0.0000000000000000
16.2500000000000000 0.0000000000000000
16.3500000000000014 0.0000000000000000
16.4499999999999993 0.0000000000000000
16.5500000000000007 0.0000000000000000
16.6500000000000021 0.0000000000000000
16.7500000000000000 0.0000000000000000
16.8500000000000014 0.0000000000000000
16.9499999999999993 0.0000000000000000
17.0500000000000007 0.0000000000000000
17.1500000000000021 0.0000000000000000
17.2500000000000000 0.0000000000000000
17.3500000000000014 0.0000000000000000
17.4499999999999993 0.0000000000000000
17.5500000000000007 0.0000000000000000
17.6500000000000021 0.0000000000000000
17.7500000000000000 0.0000000000000000
17.8500000000000014 0.0000000000000000
17.9499999999999993 0.0000000000000000
18.0500000000000007 0.0000000000000000
18.1500000000000021 0.0000000000000000
18.2500000000000000 0.0000000000000000
18.3500000000000014 0.0000000000000000
18.4499999999999993 0.0000000000000000
18.5500000000000007 0.0000000000000000
18.6500000000000021 0.0000000000000000
18.7500000000000000 0.0000000000000000
18.8500000000000014 0.0000000000000000
18.9499999999999993 0.0000000000000000
19.0500000000000007 0.0000000000000000
19.1500000000000021 0.0000000000000000
19.2500000000000000 0.0000000000000000
19.3500000000000014 0.0000000000000000
19.4500000000000028 0.0000000000000000
19.5500000000000007 0.0000000000000000
19.6500000000000021 0.0000000000000000
19.7500000000000000 0.0000000000000000
19.8500000000000014 0.0000000000000000
19.9500000000000028 0.0000000000000000
20.0500000000000007 -0.0000000000000000

dibas · Sep 26, 2012

Hi mikrom,
I have tried the idea you suggested above.
But your post has highlighted the mistake
I made of not using the appropriate
indexing of the input data and the point of
interest. I'll follow your suggestions and try
to make my procedure work.

mikrom · Sep 26, 2012

In the meantime I tried this:
lagrange.95

Code:

[COLOR=#a020f0]program[/color] Lagrange_int
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#0000ff]! A program that approximates a value of a function/experimental quantity [/color]
  [COLOR=#0000ff]! at a point based on a given set of data[/color]

  [COLOR=#0000ff]! number of data points:length of each data vector[/color]
  [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]31[/color]
  [COLOR=#0000ff]! vectors of known data: x,y-values[/color]
  [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n) :: x, y
  [COLOR=#2e8b57][b]integer[/b][/color] :: i, pos, pn
  [COLOR=#2e8b57][b]double precision[/b][/color] :: z, h, L

  [COLOR=#0000ff]! Order of Interpolation Polynomial[/color]
  pn [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]3[/color]
  
  [COLOR=#0000ff]! The data[/color]
  x [COLOR=#804040][b]=[/b][/color] ([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0.00[/color], [COLOR=#ff00ff]0.10[/color], [COLOR=#ff00ff]0.20[/color], [COLOR=#ff00ff]0.30[/color], [COLOR=#ff00ff]0.40[/color], [COLOR=#ff00ff]0.50[/color], [COLOR=#ff00ff]0.60[/color], [COLOR=#ff00ff]0.70[/color], [COLOR=#ff00ff]0.80[/color], [COLOR=#ff00ff]0.90[/color], [COLOR=#804040][b]&[/b][/color]
        [COLOR=#ff00ff]1.00[/color], [COLOR=#ff00ff]1.10[/color], [COLOR=#ff00ff]1.20[/color], [COLOR=#ff00ff]1.30[/color], [COLOR=#ff00ff]1.40[/color], [COLOR=#ff00ff]1.50[/color], [COLOR=#ff00ff]1.60[/color], [COLOR=#ff00ff]1.70[/color], [COLOR=#ff00ff]1.80[/color], [COLOR=#ff00ff]1.90[/color], [COLOR=#804040][b]&[/b][/color]
        [COLOR=#ff00ff]2.00[/color], [COLOR=#ff00ff]2.10[/color], [COLOR=#ff00ff]2.20[/color], [COLOR=#ff00ff]2.30[/color], [COLOR=#ff00ff]2.40[/color], [COLOR=#ff00ff]2.50[/color], [COLOR=#ff00ff]2.60[/color], [COLOR=#ff00ff]2.70[/color], [COLOR=#ff00ff]2.80[/color], [COLOR=#ff00ff]2.90[/color], [COLOR=#ff00ff]3.00[/color][COLOR=#804040][b]/[/b][/color])

  y [COLOR=#804040][b]=[/b][/color] ([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0.00[/color], [COLOR=#ff00ff]0.29[/color], [COLOR=#ff00ff]0.56[/color], [COLOR=#ff00ff]0.81[/color], [COLOR=#ff00ff]1.04[/color], [COLOR=#ff00ff]1.25[/color], [COLOR=#ff00ff]1.44[/color], [COLOR=#ff00ff]1.61[/color], [COLOR=#ff00ff]1.76[/color], [COLOR=#ff00ff]1.89[/color], [COLOR=#804040][b]&[/b][/color]
        [COLOR=#ff00ff]2.00[/color], [COLOR=#ff00ff]2.09[/color], [COLOR=#ff00ff]2.16[/color], [COLOR=#ff00ff]2.21[/color], [COLOR=#ff00ff]2.24[/color], [COLOR=#ff00ff]2.25[/color], [COLOR=#ff00ff]2.24[/color], [COLOR=#ff00ff]2.21[/color], [COLOR=#ff00ff]2.16[/color], [COLOR=#ff00ff]2.09[/color], [COLOR=#804040][b]&[/b][/color]
        [COLOR=#ff00ff]2.00[/color], [COLOR=#ff00ff]1.89[/color], [COLOR=#ff00ff]1.76[/color], [COLOR=#ff00ff]1.61[/color], [COLOR=#ff00ff]1.44[/color], [COLOR=#ff00ff]1.25[/color], [COLOR=#ff00ff]1.04[/color], [COLOR=#ff00ff]0.81[/color], [COLOR=#ff00ff]0.56[/color], [COLOR=#ff00ff]0.29[/color], [COLOR=#ff00ff]0.00[/color][COLOR=#804040][b]/[/b][/color])
  
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color], [COLOR=#ff00ff]'(a3,2a10)'[/color]) [COLOR=#ff00ff]'i'[/color], [COLOR=#ff00ff]'x'[/color], [COLOR=#ff00ff]'y=f(x)'[/color]
  [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color], [COLOR=#ff00ff]'(i3, f10.2 f10.2)'[/color]) i, x(i), y(i)
  [COLOR=#804040][b]enddo[/b][/color]

  [COLOR=#0000ff]! the point of evaluation[/color]
  z [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2.71[/color]

  [COLOR=#0000ff]! find the position of z in array x[/color]
  pos [COLOR=#804040][b]=[/b][/color] find_idx(x, z)

  [COLOR=#0000ff]! call function Lagrange[/color]
  [COLOR=#804040][b]if[/b][/color] (pos [COLOR=#804040][b]==[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]100[/color]) [COLOR=#804040][b]then[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Interpolation is not possible, data is out of range !'[/color]
    [COLOR=#804040][b]goto[/b][/color] [COLOR=#ff00ff]99[/color] 
  [COLOR=#804040][b]else[/b][/color] [COLOR=#804040][b]if[/b][/color] (pos [COLOR=#804040][b]==[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]then[/b][/color]
    L [COLOR=#804040][b]=[/b][/color] Lagrange(pn, x([COLOR=#ff00ff]1[/color]:pn[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color]), y([COLOR=#ff00ff]1[/color]:pn[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color]), z)
  [COLOR=#804040][b]else[/b][/color] [COLOR=#804040][b]if[/b][/color] (pos [COLOR=#804040][b]==[/b][/color] n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]then[/b][/color]
    L [COLOR=#804040][b]=[/b][/color] Lagrange(pn, x(n[COLOR=#804040][b]-[/b][/color]pn:n), y(n[COLOR=#804040][b]-[/b][/color]pn:n), z)
  [COLOR=#804040][b]else[/b][/color]
    L [COLOR=#804040][b]=[/b][/color] Lagrange(pn, x(pos[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]:pos[COLOR=#804040][b]+[/b][/color]pn[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]), y(pos[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]:pos[COLOR=#804040][b]+[/b][/color]pn[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]), z)
  [COLOR=#804040][b]end if[/b][/color]

  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]10[/color]) z, L
  [COLOR=#6a5acd]10[/color] [COLOR=#804040][b]format[/b][/color]([COLOR=#ff00ff]'The approximate value of y('[/color],[COLOR=#008080]f10.4[/color],[COLOR=#ff00ff]')='[/color],[COLOR=#008080]f10.4[/color])

[COLOR=#804040][b]  99 continue[/b][/color]
  [COLOR=#0000ff]! end  [/color]

[COLOR=#a020f0]contains[/color]
  [COLOR=#2e8b57][b]double precision[/b][/color] [COLOR=#a020f0]function[/color] Lagrange(pn, xvect, yvect, z)
    [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
    
    [COLOR=#0000ff]! Order of Interpolation Polynomial[/color]
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: pn
    [COLOR=#0000ff]! vectors of known data: x,y-values[/color]
    [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: xvect, yvect
    [COLOR=#0000ff]! the set evaluation point "z" [/color]
    [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: z

    [COLOR=#2e8b57][b]integer[/b][/color] :: i, j, n
    [COLOR=#0000ff]! Lagrange cardinal function[/color]
    [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](pn[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color]) :: L
    [COLOR=#0000ff]! the value of the polynomial/function/quantity at that point[/color]
    [COLOR=#2e8b57][b]double precision[/b][/color] :: Pz

    [COLOR=#0000ff]! n = number of data points:length of each data vector[/color]
    n [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]size[/color](xvect)
    
    [COLOR=#804040][b]if[/b][/color] (n [COLOR=#804040][b]/=[/b][/color] [COLOR=#008080]size[/color](yvect)) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Input Error !'[/color]
      [COLOR=#804040][b]return[/b][/color]
    [COLOR=#804040][b]end if[/b][/color]

    [COLOR=#0000ff]! Initializations of Pz and L[/color]
    Pz [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0.0d0[/color] [COLOR=#0000ff]! initializing the polynomia value at z[/color]
    L  [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1.0d0[/color]  [COLOR=#0000ff]! initalizing the vector of cardinal functions to 1[/color]

    [COLOR=#0000ff]! Performing the interpolation[/color]
    [COLOR=#804040][b]do[/b][/color] i [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color], n
      [COLOR=#804040][b]do[/b][/color] j [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color], n
        [COLOR=#804040][b]if[/b][/color] (i [COLOR=#804040][b]/=[/b][/color] j) [COLOR=#804040][b]then[/b][/color]
          [COLOR=#0000ff]! part of L(i)[/color]
          L(i) [COLOR=#804040][b]=[/b][/color] ( (z [COLOR=#804040][b]-[/b][/color] xvect(j)) [COLOR=#804040][b]/[/b][/color] (xvect(i) [COLOR=#804040][b]-[/b][/color] xvect(j)) )[COLOR=#804040][b]*[/b][/color] L(i)
        [COLOR=#804040][b]end if[/b][/color]
      [COLOR=#804040][b]end do[/b][/color]
      Pz [COLOR=#804040][b]=[/b][/color] Pz [COLOR=#804040][b]+[/b][/color] L(i)[COLOR=#804040][b]*[/b][/color]yvect(i) [COLOR=#0000ff]! update Pz ~ f(z)[/color]
    [COLOR=#804040][b]end do[/b][/color]

    [COLOR=#0000ff]! return value[/color]
    Lagrange [COLOR=#804040][b]=[/b][/color] Pz
  [COLOR=#a020f0]end function[/color] Lagrange

  [COLOR=#2e8b57][b]integer[/b][/color] [COLOR=#a020f0]function[/color] find_idx(xvect, x)
    [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
    [COLOR=#0000ff]! vector of x-values[/color]
    [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: xvect
    [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: x

    [COLOR=#2e8b57][b]integer[/b][/color] :: i, k, n

    [COLOR=#0000ff]! array dimension[/color]
    n [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]size[/color](xvect)

    k [COLOR=#804040][b]=[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]100[/color] [COLOR=#0000ff]! k==-100 signals ERROR[/color]
    [COLOR=#0000ff]! find index k so that x[k] < x < x[k+1][/color]
    [COLOR=#804040][b]do[/b][/color] i [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
      [COLOR=#804040][b]if[/b][/color] ((xvect(i) [COLOR=#804040][b]<=[/b][/color] x) [COLOR=#804040][b].and.[/b][/color] (x [COLOR=#804040][b]<=[/b][/color] xvect(i[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color]))) [COLOR=#804040][b]then[/b][/color]
        k [COLOR=#804040][b]=[/b][/color] i
        [COLOR=#804040][b]exit[/b][/color]
      [COLOR=#804040][b]end if[/b][/color]
    [COLOR=#804040][b]end do[/b][/color]

    find_idx [COLOR=#804040][b]=[/b][/color] k
  [COLOR=#a020f0]end function[/color] find_idx
[COLOR=#a020f0]end program[/color] Lagrange_int

Output:

Code:

$ lagrange
  i         x    y=f(x)
  1      0.00      0.00
  2      0.10      0.29
  3      0.20      0.56
  4      0.30      0.81
  5      0.40      1.04
  6      0.50      1.25
  7      0.60      1.44
  8      0.70      1.61
  9      0.80      1.76
 10      0.90      1.89
 11      1.00      2.00
 12      1.10      2.09
 13      1.20      2.16
 14      1.30      2.21
 15      1.40      2.24
 16      1.50      2.25
 17      1.60      2.24
 18      1.70      2.21
 19      1.80      2.16
 20      1.90      2.09
 21      2.00      2.00
 22      2.10      1.89
 23      2.20      1.76
 24      2.30      1.61
 25      2.40      1.44
 26      2.50      1.25
 27      2.60      1.04
 28      2.70      0.81
 29      2.80      0.56
 30      2.90      0.29
 31      3.00      0.00
The approximate value of y(    2.7100)=    0.7859

You can modify the program, so it reads the x-vect and y-vect from one file,
the z-vect from other file and for every given z-value computes y(z) using the interpolation procedure.

dibas · Sep 27, 2012

Hi mikrom,
Thanks for the modifications and the extra index-finding functions you introduced.
I did the adaptations to my problem and everything works perfectly well! Thank you
so much again.

mikrom · Sep 27, 2012

Hi dibas,
I'm glad, I could help you

Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Lagrange interpolation over a larger set of data 1

dibas

Programmer

gummibaer

Programmer

mikrom

Programmer

mikrom

Programmer

dibas

Programmer

dibas

Programmer

mikrom

Programmer

dibas

Programmer

mikrom

Programmer

Similar threads

Part and Inventory Search

Sponsor