'******************************************************************
'*   Program to demonstrate Lagrange derivative interpolation     *
'* -------------------------------------------------------------- *
'*     Ref.: Basic Scientific Subroutines Vol. II page 314        *
'*     By F.R. Ruckdeschel. BYTE/McGRAW-HILL, 1981 [1].           *
'*                                                                *
'* -------------------------------------------------------------- *
'* SAMPLE RUN:                                                    *
'*                                                                *
'*   X    2COS(2X)     YP         YP1        ERROR 1     ERROR 2  *
'* -------------------------------------------------------------- *
'* 0.00   2.000000   1.999961   2.006643  -0.0000394   0.0066434  *
'* 0.05   1.990008   1.989970   1.996569  -0.0000385   0.0065603  *
'* 0.10   1.960133   1.960096   1.966545  -0.0000373   0.0064118  *
'* 0.15   1.910673   1.910637   1.916872  -0.0000357   0.0061991  *
'* 0.20   1.842122   1.842088   1.848047  -0.0000337   0.0059246  *
'* 0.25   1.755165   1.755134   1.760756  -0.0000314   0.0055908  *
'* 0.30   1.650671   1.650642   1.655872  -0.0000288   0.0052011  *
'* 0.35   1.529684   1.529659   1.534444  -0.0000259   0.0047595  *
'* 0.40   1.393413   1.393391   1.397684  -0.0000227   0.0042704  *
'* 0.45   1.243220   1.243201   1.246958  -0.0000193   0.0037386  *
'* 0.50   1.080605   1.080589   1.083774  -0.0000157   0.0031694  *
'* 0.55   0.907192   0.907180   0.909761  -0.0000120   0.0025685  *
'* 0.60   0.724715   0.724707   0.726658  -0.0000081   0.0019420  *
'* 0.65   0.534998   0.534993   0.536294  -0.0000042   0.0012961  *
'* 0.70   0.339934   0.339934   0.340571  -0.0000002   0.0006372  *
'* 0.75   0.141474   0.141478   0.141446   0.0000038  -0.0000280  *
'* 0.80  -0.058399  -0.058391  -0.059092   0.0000078  -0.0006929  *
'* 0.85  -0.257689  -0.257677  -0.259040   0.0000116  -0.0013510  *
'* 0.90  -0.454404  -0.454389  -0.456400   0.0000154  -0.0019955  *
'* 0.95  -0.646579  -0.646560  -0.649199   0.0000190  -0.0026201  *
'* 1.00  -0.832294  -0.832271  -0.835512   0.0000224  -0.0032185  *
'* -------------------------------------------------------------- *
'******************************************************************
defint i-n
defdbl a-h,o-z

dim X(30),Y(30),XL(10),XM(10,10)

  'Derivative of sin(x) between 0 and 1
  cls
  n=4       'level of Lagrange interpolation
  ndata=26  'number of table points
  pas=0.05  'step for x
  'building X & Y Tables
  xx=0
  for i=1 to ndata
    X(i)=xx
    Y(i)=SIN(xx)
    xx = xx + pas
  next i
  'write header
  print "  X      COS(X)      YP         YP1       ERROR 1      ERROR 2  "  'heading
  print "----------------------------------------------------------------"
  'define display format
  F$="#.##  ##.######  ##.######  ##.######  ##.########  ##.########"
  'main loop of derivation
  for xx=0.05 to 1.05 step pas
    gosub 1000  'Lagrange derivative  order=4
    gosub 2000  'Parabolic derivative order=2
    print using F$; xx; COS(xx); Yp; Yp1; Yp-COS(xx); Yp1-COS(xx)
  next xx
  'write footer
  print "----------------------------------------------------------------"
  'exit program
end


'*****************************************************
'* Lagrange derivative interpolation procedure Deriv *
'* NL is the level of the interpolation ( ex. NL=3 ) *
'* N is the total number of table values.            *
'* X[i], Y[i] are the coordinate table values, Y[i]  *
'* being the dependant variable. The X[i] may be     *
'* arbitrarily spaced. x1 is the interpolation point *
'* which is assumed to be in the interval with at    *
'* least one table value to the left, and NL to the  *
'* right. Yp is returned as the desired derivative.  *
'*****************************************************
 1000 x1=xx : nl = n
    'x1 not in interval (1:ndata-nl)
    if x1>X(1) then goto 1100
    print " Error: x lower then X(1)"
    return
1100 if x1<X(ndata-nl) then goto 1200
    print " Error: x >= X(n-4)"
    return
1200 i=0
1300 i=i+1
    if x1>X(i) then goto 1300
    i=i-1
    for j=0 to nl
      XL(j)=0
      for k=0 to nl
        XM(j,k)=1
      next k
    next j
    Yp=0
    for k=0 to nl
      for j=0 to nl
        if j<>k then
          for ll=0 to nl
            if ll<>k then
              if ll=j then
                XM(ll,k)=XM(ll,k)/(X(i+k)-X(i+j))
              else
                XM(ll,k)=XM(ll,k)*(x1-X(j+i))/(X(i+k)-X(i+j))
              end if
            end if
          next ll
        end if
      next j
      for ll=0 to nl
        if ll<>k then
          XL(k)=XL(k)+XM(ll,k)
        end if
      next ll
      Yp=Yp+XL(k)*Y(i+k)
    next k
    return
  end


'*************************************************************
'* CALCULE LES COEFFICIENTS A,B DE LA PARABOLE Y=A*X*X+B*X+C *
'* PASSANT PAR LES 3 POINTS : (X1,Y1),(X2,Y2) ET (X3,Y3)     *
'* COEFFICIENT C NON UTILISE ICI.                            *
'* --------------------------------------------------------- *
'* Calculates coefficients, a, b of parabola Y=A*X+X+B*X+C   *
'* passing through 3 points: (X1,Y1), (X2,Y2) and (X3,Y3).   *
'* Coefficient c is not used here.                           *
'*************************************************************
1900 alpha=xx1*xx1-x2*x2  'x1 is called here xx1
  beta=xx1-x2
  gamma=x2*x2-x3*x3
  delta=x2-x3
  a=(y1-2#*y2+y3)/(alpha-gamma)
  b=(y1-y2-alpha*a)/beta
  return
End

'****************************************************
'*     Interpolation of order=2 ( parabola )        *
'*                            by J-P Moreau         *
'****************************************************
2000 x1 = xx
  if x1>X(1) then goto 2100
  print " Error in order 2" : return
2100 if x1<X(ndata-2) then goto 2200
  print " Error in order 2" : return
2200 i=0
2300 i=i+1
  if x1>X(i) then goto 2300
  i=i-1
  xx1=X(i):y1=Y(i):x2=X(i+1):y2=Y(i+1):x3=X(i+2):y3=Y(i+2)
  gosub 1900
  'estimation of derivative in x1
  Yp1=2#*a*x1+b
  return
End

' End of file derivati.bas