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:

  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes" - Cambridge University Press