Euler-Lagrange Equation for the HP-41

Overview
 

 1°)  First Order Lagrangian - One Function

  a)  Without Calculus
  b)  With Calculus

 2°)  Second Order Lagrangian - One Function

  a)  Without Calculus
  b)  With Calculus

 3°)  First Order Lagrangian - Two Functions

  a)  Without Calculus
  b)  With Calculus
 

-These programs solve problems of variational calculus of the type: find an extremum of  the integral     I = §_a^b  L dx

  where L is a function of x , y , y' in paragraph 1 , a function of  x , y , y' , y"  in paragraph 2  and a function of  x , y , z , y' , z'  in paragraph 3

-In paragraphs 1a) 2a) 3a) , you only have to program the function L
-In paragraphs 1b) 2b) 3b) , you also have to calculate and program the Euler-Lagrange equation(s)

-In the first case, the derivatives are approximated by 2nd order formulae.
-In the second case, the Euler-Lagrange equation is solved by a Runge-Kutta method: the precision is much better and the routine is much faster !
 

1°)  First Order Lagrangian - One Function
 

     a)  Without Calculus
 

-The following programs try to find an extremum of the integral     I = §_a^b  L( x , y , y' ) dx       ( § = integral )

   where y(x) is an unknown function satisfying    y(a) = A , y(b) = B

-The Euler-Equation is   ¶L / ¶y = (d/dx) ( ¶L / ¶y' )                where  ¶ = partial derivative

   whence     ¶L / ¶y = ¶²L / ¶x ¶y' + y' ¶²L / ¶x ¶y' + y'' ¶²L / ¶y'²
 

>>> Here, we assume that  ¶²L / ¶y'²  # 0    so,    y'' =  ( ¶L / ¶y - ¶²L / ¶x ¶y' - y' ¶²L / ¶x ¶y' ) / ( ¶²L / ¶y'² )

-In this routine, the derivatives are approximated by 2nd order formulas and the integral is computed by Simpson's rule.
-The value of y'(a) is found by the secant method ( lines 19 to 40 )
 
 

Data Registers:           •  R00 = N = a positive even integer = Number of subintervals

                                       ( Registers R00 thru R04 are to be initialized before executing "EULAG" )

                                      •  R01 = a     •  R03 = y(a)       R05 = y'(a)      R07 =  I_min = §_a^b  L( x , y , y' ) dx     R09 to R19: temp
                                      •  R02 = b     •  R04 = y(b)      R06 = y'(b)      R08 =  20.eee

        And if  F00 is clear:    R20 = y(a)  R21 = y(a+h)  R22 = y(a+2.h)  .......................   Reee = y(b)       with  h = (b-a) / N

Flags:  F00 & F10       CF 00  ->  the values of the function y(x)  are stored in  R20 , R21 , .....
                                     SF 00  ->  no !

Subroutine:   A program, called  LBL "L"  calculating  L( x , y , y' )  in X-register with x , y , y'  in registers  X , Y , Z  respectively.
 

-Line 207 is a three-byte GTO
 
 

 
     ( 338 bytes / SIZE 021 or 022 + N )
 
 

   Where  m & n are your initial guesses for y'(a) ,  I = §_a^b  L( x , y , y' ) dx

Example:     Minimize  I  =  §₁²  ( x² + y² + y'² ) dx  with  y(1) = y(2) = 1

-Load in main memory:
 
 

 
-Store a , b , y(a) , y(b)  into  R01 thru R04:

   1  STO 01    1  STO 03  STO 04
   2  STO 02

-If you choose  N = 10    and the initial guesses 1 & -1 for y'(a)

   10  STO 00   FIX 4

    1   ENTER^
   CHS  XEQ "EULAG"  >>>>      I_min   =  3.2575                                            ---Execution time = 5mn33s---
                                       RDN      y'(1)   = -0.4661
                                       RDN      y'(2)   =  0.4587
                                       RDN    20.eee  = 20.0300
 

-And we find in R20 to R30:
 
 

 
-With N = 50 & FIX 6  it yields  I = 3.257563
-With N = 100 & FIX 9 ------   I = 3.257566529

-The exact solution is y(x) = [ Cosh ( x - 3/2 ) ] / Cosh (1/2)
 
 

 
-And the exact minimum integral is  I_min = 3.257567648  ( rounded to 9 decimals )

Notes:

-The precision is controlled by the display format ( line 38 )
-The HP41 displays the successive approximations of y'(a) in R05

-However, the derivatives are not calculated by very accurate formulae and the differential equation is solved by

    y(x+h) ~ y(x) + h.y'(x) + (h²/2).y"(x)  &  y'(x+h) ~ y'(x) + h.y"(x)

 so we cannot expect a great precision...
 

     b)  With Calculus
 

-Before using this program, you have to compute  y" = T(x,y,y')
-Thus, we can use a Runge-Kutta formula to find the solution, with a much better precision.
 

Data Registers:           •  R00 = N = a positive even integer = Number of subintervals

                                       ( Registers R00 thru R04 are to be initialized before executing "EULAG" )

                                      •  R01 = a     •  R03 = y(a)       R05 = y'(a)      R07 =  I_min = §_a^b  L( x , y , y' ) dx     R09 to R19: temp
                                      •  R02 = b     •  R04 = y(b)      R06 = y'(b)      R08 =  20.eee

        And if  F00 is clear:    R20 = y(a)  R21 = y(a+h)  R22 = y(a+2.h)  .......................   Reee = y(b)       with  h = (b-a) / N

Flags:  F00 & F10       CF 00  ->  the values of the function y(x)  are stored in  R20 , R21 , .....
                                     SF 00  ->  no !

Subroutines:   A program, called  LBL "L"  calculating  L( x , y , y' )  in X-register with x , y , y'  in registers  X , Y , Z  respectively.
                         and another one     LBL "T"  ----------  y" = T(x,y,y') --------------------------------------------------------------
 

-Line 147 is a three-byte GTO 03
 
 

 
     ( 254 bytes / SIZE 021 or 022 + N )
 
 

   Where  m & n are your initial guesses for y'(a) ,  I = §_a^b  L( x , y , y' ) dx

Example:     Minimize  I  =  §₁²  ( x² + y² + y'² ) dx  with  y(1) = y(2) = 1

-Here, Euler-Lagrange equation gives:   ¶L / ¶y = (d/dx) ( ¶L / ¶y' )  --->  2.y = (d/dx) (2.y')

   whence   y" = y

-Load in main memory:
 
 

 
-Store a , b , y(a) , y(b)  into  R01 thru R04:

   1  STO 01    1  STO 03  STO 04
   2  STO 02

-If you choose  N = 10    and the initial guesses 1 & -1 for y'(a)

   10  STO 00   FIX 6

    1   ENTER^
   CHS  XEQ "EULAG"  >>>>      I_min   =  3.257576                                            ---Execution time = 1mn58s---
                                       RDN      y'(1)   = -0.462117
                                       RDN      y'(2)   =  0.462117
                                       RDN    20.eee  = 20.0300
 

-And we find in R20 to R30:
 
 

 
-With N = 50 & FIX 9  it yields  I = 3.257567662
-With N = 100 & FIX 9 ------   I = 3.257567651

-The exact solution is y(x) = [ Cosh ( x - 3/2 ) ] / Cosh (1/2)
-So the errors are smaller than E-6
-The exact values of y'(a) & y'(b) are -/+ 0.462117157...

-And the exact minimum integral is  I_min = 3.257567648  ( rounded to 9 decimals )

Notes:

-The precision is controlled by the display format ( line 39 )
-The HP41 displays the successive approximations of y'(a) in R05
 

2°)  Second Order Lagrangian - One Function
 

     a)  Without Calculus
 

-Here, we try to find the extremum of  I = §_a^b  L( x , y , y' , y" ) dx   with  y(a) = A , y(b) = B , y'(a) = A' , y'(b) = B'

-In this case, the Euler-Lagrange equation becomes:   ¶L / ¶y - (d/dx) ( ¶L / ¶y' ) + (d²/dx²) ( ¶L / ¶y'' ) = 0            where  ¶ = partial derivative
 

>>> We assume that  ¶²L / ¶y''²  # 0
 

-The derivatives are also approximated by second order formulae.
-The values of y"(a) & y'''(a) are found by a 2-dimensional secalt method ( lines 20 to 90 )
 
 

Data Registers:           •  R00 = N = a positive even integer = Number of subintervals

                                       ( Registers R00 thru R06 are to be initialized before executing "2ELG" )

                                      •  R01 = a     •  R03 = y(a)     •  R05 = y'(a)     R07 = y"(a)   R09 = y'''(a)   R11 =  I_min = §_a^b  L( x , y , y' , y" ) dx
                                      •  R02 = b     •  R04 = y(b)     •  R06 = y'(b)    R08 = y"(b)   R10 = y'''(b)   R12 = 33.eee    R13 to R32: temp

        And if  F00 is clear:    R33 = y(a)  R34 = y(a+h)  R35 = y(a+2.h)  .......................   Reee = y(b)       with  h = (b-a) / N

Flags:  F00 & F02-F03-F10       CF 00  ->  the values of the function y(x)  are stored in  R33 , R34 , .....
                                                    SF 00  ->  no !

Subroutine:   A program, called  LBL "L"  calculating  L( x , y , y' ,y" )  in X-register with x , y , y' , y"  in registers  X , Y , Z , T  respectively.
 
 

-Lines 100 & 630 are three-byte GTOs
 
 

 
        ( 985 bytes SIZE 033 or 34+N )
 
 

   Where  m , n & m' , n' are your initial guesses for y''(a) & y'''(a)  and  I = §_a^b  L( x , y , y' , y" ) dx

Example:     Minimize  I  =  §₀¹  ( x + y² + 2 y'² + y''² ) dx  with  y(0) = 0 ,  y(1) = Cosh 1 , y'(0) = 1 , y'(1) = e

-Load in main memory:
 
 

 
-Store a , b , y(a) , y(b) , y'(a) , y'(b)  into  R01 thru R06:

   0  STO 01              0           STO 03    1           STO 05
   1  STO 02    1.543080635  STO 04    1  E^X  STO 06

-If you choose  N = 10    and the initial guesses 0 , 0.1 for y''(0)  and the same for y'''(0)

   10  STO 00   FIX 4

    0       ENTER^
   0.1     ENTER^
    0       ENTER^
   0.1     XEQ "2ELG"  >>>>      I_min   =  10.5164                                        ---Execution time = 8s---               with V41 in turbo mode
                                    RDN      y"(0)   = -0.0230 = R07
                                    RDN      y"(1)   =   3.8681 = R08
                                    RDN    33.eee  = 33.043

-In R09 & R10 we have  y'''(0) = 3.1340  &  y'''(1) = 5.6297

-And we find in R33 to R43:
 
 

 
-With N =   50   & FIX 6 it yields  I = 10.515922
-With N = 1000 & FIX 9 -------   I = 10.515916337     ( these 2 results obtained with free42 )

-The exact solution is y(x) = x Cosh x
 
 

 
-So, y"(0) = 0 & y"(1) = 2 Sinh 1 + Cosh 1 = 3.893483022...  , y'''(0) = 3 & y'''(1) = 3 Cosh 1 + Sinh 1 = 5.804443098...

-And the exact minimum integral is  I_min = ( 3 e² - 1 - e^-2 ) / 2 = 10.515916507   ( rounded to 9 decimals )

Notes:

-The precision is controlled by the display format ( line 87 )
-The HP41 displays the successive approximations of y"(a) in R07

-This program should be used with V41 or free42 !
 

     b)  With Calculus
 

-We try to find the extremum of  I = §_a^b  L( x , y , y' , y" ) dx   with  y(a) = A , y(b) = B , y'(a) = A' , y'(b) = B'

-The Euler-Lagrange equation is:     ¶L / ¶y - (d/dx) ( ¶L / ¶y' ) + (d²/dx²) ( ¶L / ¶y'' ) = 0            where  ¶ = partial derivative
 

>>> We assume that  ¶²L / ¶y''²  # 0
 

-You also have to calculate y''''(x) as a function of  x , y , y' , y" , y'''  and write a program - named LBL "T" - to calculate y'''' = T(x,y,y',y",y''')

-This 4th-order differential equation is solved by the "classical" Runge-Kutta method.
 

Data Registers:           •  R00 = N = a positive even integer = Number of subintervals

                                       ( Registers R00 thru R06 are to be initialized before executing "2ELG" )

                                      •  R01 = a     •  R03 = y(a)     •  R05 = y'(a)     R07 = y"(a)   R09 = y'''(a)   R11 =  I_min = §_a^b  L( x , y , y' , y" ) dx
                                      •  R02 = b     •  R04 = y(b)     •  R06 = y'(b)    R08 = y"(b)   R10 = y'''(b)   R12 = 37.eee    R13 to R36: temp

        And if  F00 is clear:    R37 = y(a)  R38 = y(a+h)  R39 = y(a+2.h)  .......................   Reee = y(b)       with  h = (b-a) / N

Flags:  F00 & F10            CF 00  ->  the values of the function y(x)  are stored in  R37 , R38 , .....
                                          SF 00  ->  no !

Subroutines:   A program, called  LBL "L"  calculating  L( x , y , y' ,y" )  in X-register with x , y , y' , y"  in registers R11-R12-R13-R14  respectively.
                        and another, called LBL "T"  calculating  T( x , y , y' , y" , y''' ) in X-register with x , y , y' , y" , y''' in R11-R12-R13-R14-R15 respectively.
 
 

-Line 266 is a three-byte GTO 03
 
 

 
        ( 460 bytes SIZE 037 or 38+N )
 
 

   Where  m , n & m' , n' are your initial guesses for y''(a) & y'''(a)  and  I = §_a^b  L( x , y , y' , y" ) dx

Example:     Minimize  I  =  §₀¹  ( x + y² + 2 y'² + y''² ) dx  with  y(0) = 0 ,  y(1) = Cosh 1 , y'(0) = 1 , y'(1) = e

-Writing    ¶L / ¶y - (d/dx) ( ¶L / ¶y' ) + (d²/dx²) ( ¶L / ¶y'' ) = 0     gives:

   2.y - 4.y" + 2.y'''' = 0   whence  y'''' = 2.y" - y
 

-Load in main memory:
 
 

 
-Store a , b , y(a) , y(b) , y'(a) , y'(b)  into  R01 thru R06:

   0  STO 01              0           STO 03    1           STO 05
   1  STO 02    1.543080635  STO 04    1  E^X  STO 06

-If you choose  N = 10    and the initial guesses 0 , 0.1 for y''(0)  and the same for y'''(0)

   10  STO 00   FIX 6

    0       ENTER^
   0.1     ENTER^
    0       ENTER^
   0.1     XEQ "2ELG"  >>>>      I_min   =  10.516312                                            ---Execution time = 6m15s---
                                    RDN      y"(0)  =  0.000030 = R07
                                    RDN      y'''(0) =  2.999942 = R09
                                    RDN    37.eee  = 37.047

-In R08 & R10 we have  y"(1)   =  3.893458  &  y'''(1) = 5.804386

-And we find in R37 to R47:
 
 

 
-The exact solution is y(x) = x Cosh x
 
 

 
-Of course, the precision is much better !
 

-With N =  40  & FIX 6 , it yields:   I = 10.515918
-With N = 100 & FIX 8 , it yields:   I = 10.51591654

-The exact minimum integral is  I_min = ( 3 e² - 1 - e^-2 ) / 2 = 10.515916507   ( rounded to 9 decimals )
 

Notes:

-The precision is controlled by the display format ( line 88 )
-The HP41 displays the successive approximations of y"(a) in R07
 

3°)  First Order Lagrangian - Two Functions
 

     a)  Without Calculus
 

-The Euler-Lagrange equation becomes a system of 2 differential equations:

      ¶L / ¶y = (d/dx) ( ¶L / ¶y' )
      ¶L / ¶z = (d/dx) ( ¶L / ¶z' )

-After some calculus, it yields:

     y" ¶²L / ¶y'² + z" ¶²L / ¶y'¶z' = ¶L / ¶y - ¶²L / ¶x¶y' - y' ¶²L / ¶y¶y' - z' ¶²L / ¶z¶y'

     y" ¶²L / ¶y'¶z' + z" ¶²L / ¶z'² = ¶L / ¶z - ¶²L / ¶x¶z' - y' ¶²L / ¶y¶z' - z' ¶²L / ¶z¶z'
 

>>>> Assuming the determinant    ( ¶²L / ¶y'² ) ( ¶²L / ¶z'² ) - ( ¶²L / ¶y'¶z' )²  #  0     the following routine solves this system.
 

Data Registers:           •  R00 = N = a positive even integer = Number of subintervals

                                       ( Registers R00 thru R06 are to be initialized before executing "ELG2" )

                                      •  R01 = a     •  R03 = y(a)     •  R05 = z(a)     R07 = y'(a)   R09 = z'(a)   R11 =  I_min = §_a^b  L( x , y , z , y' , z' ) dx
                                      •  R02 = b     •  R04 = y(b)     •  R06 = z(b)    R08 = y'(b)   R10 = z'(b)   R12 = 43.eee    R13 to R42: temp

        And if  F00 is clear:    R43 = y(a)  R44 = y(a+h)  R45 = y(a+2.h)  .......................   Reee = y(b)       with  h = (b-a) / N
                                          Ree+1 = z(a)  .............................................................  Ree+N+1 = z(b)

Flags:  F00 & F10            CF 00  ->  the values of the function y(x)  are stored in  R43 , R44 , .....
                                          SF 00  ->  no !

Subroutine: A program, called  LBL "L"  calculating  L( x , y , z , y' , z' )  in X-register with x , y , z , y' , z'  in registers R11-R12-R13-R14-R15  respectively.
 

-Line 325 is a three-byte GTO 03
 
 

 
        ( 605 bytes SIZE 045 or 46+2.N )
 
 

   Where  m , n & m' , n' are your initial guesses for y'(a) & z'(a)  and  I = §_a^b  L( x , y , z , y' , z' ) dx

Example:     Minimize  I  =  §₁²  ( x y² z - x² z² + y'² + z'² / 2 ) dx   with  y(1) = 1 ,  y(2) = 2 ,  z(1) = 0 ,  z(2) = 1

-Load this subroutine in main memory:
 
 

 
-Store a , b , y(a) , y(b) , z(a) , z(b)  into  R01 thru R06:

   1  STO 01   1   STO 03   0    STO 05
   2  STO 02   2   STO 04   1    STO 06

-With  N = 10   10  STO 00  and if you choose 0  1  as initial guesses for  y'(1) & z'(1)

        FIX 4  CF 00

    0  ENTER^
    1  ENTER^
    0  ENTER^
    1  XEQ "ELG2"  >>>>       I_min   =  2.7454                                           ---Execution time = 9s---                     With V41 in turbo mode
                               RDN      y'(1)   =  0.7864 = R07
                               RDN      z'(1)   =  0.3864 = R09
                               RDN    43.eee  = 43.0530

-In R08 & R10 we have  y'(2) = 1.7293  &  z'(2) = 1.4937

-And we find in R43 to R53 & R54 to R64
 
 

 
-With  N = 100 & FIX 6   V41 returns  I = 2.746811  ( free42 gives  2.746812734 )
-With  N = 1000 , free42 ->  2.746832062

-Compared to the results obtained by the following program, these values are not very accurate...
 

Notes:

-The precision is controlled by the display format ( line 99 )
-The HP41 displays the successive approximations of y'(a) in R07
 

     b)  With Calculus
 

-The Euler-Lagrange equations are:

      ¶L / ¶y = (d/dx) ( ¶L / ¶y' )
      ¶L / ¶z = (d/dx) ( ¶L / ¶z' )

-This system of 2 differential equations is solved by the "classical" Runge-Kutta method, assuming that we may express

      y" = T(x,y,z,y',z')   &  z" = U(x,y,z,y',z')

-You have to find these 2 differential equtions and write a program to compute T & U
 
 

Data Registers:           •  R00 = N = a positive even integer = Number of subintervals

                                       ( Registers R00 thru R06 are to be initialized before executing "ELG2" )

                                      •  R01 = a     •  R03 = y(a)     •  R05 = z(a)     R07 = y'(a)   R09 = z'(a)   R11 =  I_min = §_a^b  L( x , y , z , y' , z' ) dx
                                      •  R02 = b     •  R04 = y(b)     •  R06 = z(b)    R08 = y'(b)   R10 = z'(b)   R12 = 37.eee    R13 to R36: temp

        And if  F00 is clear:    R37 = y(a)  R38 = y(a+h)  R39 = y(a+2.h)  .......................   Reee = y(b)       with  h = (b-a) / N
                                          Ree+1 = z(a)  .............................................................  Ree+N+1 = z(b)

Flags:  F00 & F10            CF 00  ->  the values of the function y(x)  are stored in  R37 , R38 , .....
                                          SF 00  ->  no !

Subroutines: A program, called  LBL "L"  calculating  L( x , y , z , y' , z' )  in X-register with x , y , z , y' , z'  in registers R11-R12-R13-R14-R15  respectively.
                        and another, called LBL "T"  calculating  T( x , y , z , y' , z' ) in Y-register
                                                                             and    U( x , y , z , y' , z' ) in X-register with x , y , z , y' , z'  in R11-R12-R13-R14-R15 respectively.
 
 

-Line 278 is a three-byte GTO 03
 
 

 
        ( 482 bytes SIZE 038 or 39+2.N )
 
 

   Where  m , n & m' , n' are your initial guesses for y'(a) & z'(a)  and  I = §_a^b  L( x , y , z , y' , z' ) dx

Example:     Minimize  I  =  §₁²  ( x y² z - x² z² + y'² + z'² / 2 ) dx   with  y(1) = 1 ,  y(2) = 2 ,  z(1) = 0 ,  z(2) = 1

-The Euler-Lagrange system is:

   y" =   x y z
   z" = - 2 x² z + x y²

-Load these subroutines in main memory:
 
 

 
-Store a , b , y(a) , y(b) , z(a) , z(b)  into  R01 thru R06:

   1  STO 01   1   STO 03   0    STO 05
   2  STO 02   2   STO 04   1    STO 06

-With  N = 10      10  STO 00  and if you choose 0  1  as initial guesses for  y'(1) & z'(1)

        FIX 4  CF 00

    0  ENTER^
    1  ENTER^
    0  ENTER^
    1  XEQ "ELG2"  >>>>       I_min   =  2.7473                                                           ---Execution time = 13mn54---
                               RDN      y'(1)   =  0.7200 = R07
                               RDN      z'(1)   =  0.4538 = R09
                               RDN    37.eee  = 37.0470

-In R08 & R10 we have  y'(2) = 1.8909  &  z'(2) = 1.3289

-And we find in R37 to R47 & R48 to R58
 
 

 
-With  N = 100 & FIX 9   V41 returns  I = 2.746832313  ( free42 gives  2.746832312 )
-With  N = 1000 , free42 ->  2.746832261
 

Notes:

-The precision is controlled by the display format ( line 88 )
-The HP41 displays the successive approximations of y'(a) in R07

-I don't know the exact solution of this problem !
 

Reference:

[1] - https://en.wikipedia.org/w/index.php?title=Euler–Lagrange_equation&oldid=912228555