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  c₀ , c₁ , ............ , c_n
-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  c₀ , c₁ , ............ , c_n  may be computed by the formula:

   c_j = [2/(n+1)]  SUM_k=0,1,...,n  Cos [ 180° j ( k+1/2 ) / ( n+1 ) ]  f { Cos [ 180° ( k+1/2 ) / ( n+1 ) ] }               if  j # 0
   c_j = [1/(n+1)]  SUM_k=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) ~ c₀ + c₁ T₁(x) + ............ + c_n T_n(x)    where    T_k(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 = c₀ , R19 = c₁ , ...... , R18+n = c_n

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+x² )    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 = c₀ =  0.061130486                            R24 = c₆ =  0.000001210
  R19 = c₁ = -0.038060320                            R25 = c₇ =  0.000000410
  R20 = c₂ =  0.008422922                            R26 = c₈ = -0.000000170
  R21 = c₃ = -0.001522665                            R27 = c₉ =  0.000000042
  R22 = c₄ =  0.000227407                            R28 = c₁₀ = -0.000000008
  R23 = c₅ = -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 = c₀ , •  R19 = c₁ , ...... , •  R18+n = c_n

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+x² )    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