hp41programs

Chebapprox

Chebyshev Approximation for the HP-41


Overview
 

 1°)  Chebyshev Coefficients
 2°)  Chebyshev Approximation
 

-Given a function f defined over an interval [a,b], the 1st program calculates the Chebyshev coefficients  c0 , c1 , ............ , cn
-Then, the second routine approximate the function and its derivative.
 
 

1°)  Chebyshev Coefficients
 
 

-If f is a function defined over [-1,+1] and if n is a positive integer, the Chebyshev coefficients  c0 , c1 , ............ , cn  may be computed by the formula:

   cj = [2/(n+1)]  SUMk=0,1,...,n  Cos [ 180° j ( k+1/2 ) / ( n+1 ) ]  f { Cos [ 180° ( k+1/2 ) / ( n+1 ) ] }               if  j # 0
   cj = [1/(n+1)]  SUMk=0,1,...,n  Cos [ 180° j ( k+1/2 ) / ( n+1 ) ]  f { Cos [ 180° ( k+1/2 ) / ( n+1 ) ] }               if  j = 0

-Then,  f(x) ~ c0 + c1 T1(x) + ............ + cn Tn(x)    where    Tk(x)  are the Chebyshev polynomials of the first kind

-If f is defined over  [a,b] , we make the change of variable  y = ( 2 x - a - b ) / ( b - a )
 

Data Registers:           •  R10 = function name                ( Registers R10 thru R12 are to be initialized before executing  "CHBCF" )

                                      •  R11 = a                                     R13 to R17: temp                R00 thru R09 are unused
                                      •  R12 = b

     >>> When the program stops,   R13 = control number of the coefficients and   R18 = c0 , R19 = c1 , ...... , R18+n = cn

Flags: /
Subroutine:  A program that takes x in X-register and returns f(x) in X-register without disturbing R10 , R11 , ....
 
 

 01  LBL "CHBCF"
 02  DEG
 03  STO 13          
 04  STO 14 
 05  18.017
 06  +
 07  STO 17
 08  ISG 13
 09  LBL 01
 10  RCL 13
 11  STO 15
 12  CLX
 13  STO IND 17
 14  LBL 02
 15  RCL 15          
 16  .5
 17  -
 18  PI
 19  R-D
 20  *
 21  RCL 13
 22  /
 23  STO 16
 24  COS
 25  RCL 12
 26  RCL 11
 27  -
 28  *
 29  RCL 11          
 30  RCL 12
 31  +
 32  +
 33  2
 34  /
 35  XEQ IND 10
 36  RCL 14
 37  RCL 16
 38  *
 39  COS
 40  *
 41  ST+ IND 17
 42  DSE 15
 43  GTO 02
 44  RCL 14          
 45  X#0?
 46  SIGN
 47  1
 48  +
 49  RCL 13
 50  /
 51  ST* IND 17
 52  DSE 14
 53  CLX
 54  DSE 17
 55  GTO 01
 56  RCL 13          
 57  17
 58  +
 59   E3
 60  /
 61  18
 62  +
 63  STO 13
 64  END

 
    ( 96 bytes / SIZE 19+nnn )
 
 

      STACK        INPUT      OUTPUT
           X             n         18.eee

   where n is a positive integer and 18.eee = control number of the Chebyshev coefficients ( eee = 18+n )

Example:    f(x) = 1 / ( PI+x+x2 )    and   [a,b] =  [2,5]

-Load a routine to calculate f(x) , for instance:

 LBL "T"
 ENTER^
 X^2
 +
 PI
 +
 1/X
 RTN

 "T"  ASTO 10
  2   STO 11
  5   STO 12

-And if you choose n = 10 ( 11 coefficients )

  10  XEQ "CHBCF"   >>>>   18.028                            ---Execution time = 5m05s---

-You get in registers R18 thru R28

  R18 = c0 =  0.061130486                            R24 = c6 =  0.000001210
  R19 = c1 = -0.038060320                            R25 = c7 =  0.000000410
  R20 = c2 =  0.008422922                            R26 = c8 = -0.000000170
  R21 = c3 = -0.001522665                            R27 = c9 =  0.000000042
  R22 = c4 =  0.000227407                            R28 = c10 = -0.000000008
  R23 = c5 = -0.000025750

Notes:

-"CHBCF" runs in DEG mode, so the subroutine must restore this angular mode at the end if need be.
-Alternatively, add FC? 43 after line 18 and delete line 02: it will work in DEG or RAD mode.
 

2°)  Chebyshev Approximation
 

-After executing "CHBCF" , we can use "CdT" to approximate f(x) or f'(x)
 

Data Registers:           •  R10 = function name                ( Execute "CHBCF" before executing  "CHBAP" )

                                      •  R11 = a
                                      •  R12 = b

                                      •  R13 = control number of the coefficients and   •  R18 = c0 , •  R19 = c1 , ...... , •  R18+n = cn

Flag:  CF 01 to calculate f(x)
           SF 01 to calculate f'(x)

Subroutine:  "CdT" - M-Code routine or focal program ( cf "Ephemerides and Chebyshev Polynomials" in the "Astronomy" section ).

-If you replace line 15 by  XEQ "CdT" ( focal program ) , delete lines 12-13
 
 

 01  LBL "CHBAP"
 02  ST+ X
 03  RCL 11
 04  RCL 12
 05  +
 06  -
 07  RCL 12
 08  RCL 11
 09  -
 10  /
 11  RCL 13
 12  FS? 01
 13  CHS
 14  X<>Y
 15  CdT
 16  2
 17  RCL 12
 18  RCL 11
 19  -
 20  /
 21  X<>Y
 22  FS? 01
 23  *
 24  END

 
( 38 bytes / SIZE 19+nnn )
 
 

      STACK        INPUT      OUTPUT
           X             x     f(x) or f'(x)

  where  a <= x <= b ,  f(x) if CF 01 , f'(x) if SF 01

Example:   The same one:    f(x) = 1 / ( PI+x+x2 )    and   [a,b] =  [2,5]

-After executing "CHBCF", let's evaluate f(3) & f'(3)

  3   CF 01  XEQ "CHBAP"  >>>>  f(3) ~ 0.066043252
  3   SF 01          R/S             >>>>  f'(3) ~ -0.030531990

Notes:

-In this example, the exact values - rounded to 9D - are  f(3) = 0.066043251  &  f'(3) = -0.030531977

-Choosing a larger n-value would give a better precision.

-The Chebyshev polynomials are useful to approximate f(x) if f is very complicated - like the planetary positions.
-But it's also interesting to use these programs to evaluate the derivative f'(x):
-The results are often more accurate than those given by the other numerical differentiation methods.
 
 

Reference:

[1]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes" - Cambridge University Press