hp41programs

Clausen's Integral

Clausen's Integral for the HP-41


Overview
 

-The following program computes   f(x) = - §0x Ln ( 2 sin (t/2) ) dt                     ( § = Integral )

-Integration by parts gives:   f(x) = -x Ln ( 2 sin (x/2) ) + 2 §0x/2 ( t / Tan t ) dt
-The series expansion   t / Tan t = 1 - x2/3 - x4/45 - 2 x6/945 - ..... - [ (-1)n-1 22n/(2n)! ] B2n x2n - .....................  is used  where  B2n are Bernoulli numbers.
-Bernoulli numbers B2n are calculated through  Zeta function for n > 3
 

Program Listing
 

Data Registers:    R00-R01  are used by "ZETA"

                               R02 = x    R03 thru R06: temp   (  R05 = §0x/2 ( t / Tan t ) dt )
Flags: /
Subroutine:   "ZETA"  ( cf "Riemann Zeta Function for the HP-41" )

-If x is small enough ( for example  x = 0.15 ), "ZETA" is not needed.
 
 

01  LBL "CLAUS"
02  DEG
03  STO 02
04  2
05  /
06  ENTER^
07  X^2
08  STO 04          
09  LASTX
10  *
11  STO 03
12  9
13  /
14  -
15  RCL 03
16  RCL 04		 17  *
18  STO 03 
19  225
20  /
21  -
22  RCL 03
23  RCL 04 
24  *
25  STO 03          
26  3307.5
27  /
28  CHS
29  X<>Y
30  +
31  STO 05
32  LASTX		 33  X=Y? 
34  GTO 00
35  2
36  PI
37  6
38  STO 06 
39  Y^X
40  /
41  CHS
42  ST* 03
43  FIX 9
44  LBL 01
45  RCL 06          
46  2
47  +
48  STO 06		 49  XEQ "ZETA" 
50  RCL 03
51  RCL 04 
52  *
53  PI
54  X^2
55  /
56  STO 03          
57  *
58  RCL 06
59  1
60  +
61  /
62  RCL 05
63  +
64  STO 05		 65  LASTX
66  X#Y?
67  GTO 01
68  LBL 00
69  ST+ X
70  RCL 02          
71  2
72  /
73  R-D 
74  SIN
75  ST+ X
76  LN
77  RCL 02
78  *
79  -
80  END

 
   ( 108 bytes / SIZE 007 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             /           g(x)
           X             x           f(x)

  where  f(x) = Clausen's Integral            ( 0 < x < 2 PI )
   and    g(x) = §0x/2 ( t / Tan t ) dt

Examples:

  PI  3  /   XEQ "CLAUS"  >>>>   f(PI/3) = 1.014941606                ---Execution time = 28s---
                                          X<>Y   g(PI/3) = 0.507470803

  0.15   R/S    >>>>   f(0.15) = 0.434614878       ( in 2 seconds )
                     X<>Y   g(0.15) = 0.074953114

    6      R/S    >>>>   f(6) = -0.640782664           ( in 7mn49s )
                     X<>Y   g(6) = -4.115383667
 

Notes:

-Execution time tends to infinity as x tends to 2 PI
-This routine produces DATA ERROR if x = 0 but since f(0) = 0, simply add  X=0?  RTN  after line 01

-Clausen's integral is also equal to  SUMk=1,2,......  [ Sin (kx) ] / k2 but thousands of terms would be required to return accurate results if we used this formula.
-However, Clausen's integral Cl(x) is closely related to Lobachevsky function L(x), we have:

   Cl(x) = 2 L(x/2)     at least if  x >= 0 & x <= 2.PI

-So, we can use the program "LOB" listed in "Lobachevsky Function for the HP-41"
 

Reference:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4