Euler-Lagrange

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

01 LBL "EULAG"
02 STO 06
03 X<>Y
04 STO 05
05 RCL 00
06 2
07 MOD
08 ST+ 00
09 RCL 02
10 RCL 01
11 -
12 RCL 00
13 /
14 STO 10
15 SF 10
16 RCL 06
17 XEQ 02
18 STO 11
19 LBL 01
20 VIEW 05
21 RCL 05
22 XEQ 02
23 ENTER
24 ENTER
25 X<> 11
26 -
27 X#0?
28 /
29 RCL 06
30 RCL 05
31 STO 06
32 STO T
33 -
34 *

35 +
36 STO 05
37 -
38 RND
39 X#0?
40 GTO 01
41 20
42 STO 19
43 2
44 STO 16
45 X^2
46 STO 15
47 RCL 05
48 CF 10
49 LBL 02
50 STO 09
51 RCL 01
52 STO 07
53 RCL 00
54 STO 18
55 RCL 03
56 STO 08
57 FS? 10
58 GTO 03
59 RCL 09
60 X<>Y
61 STO 20
62 RCL 07
63 XEQ "L"
64 CHS
65 STO 14
66 LBL 03
67 RCL 09
68 RCL 08

69 RCL 10
70 +
71 RCL 07
72 XEQ "L"
73 STO 12
74 RCL 09
75 RCL 08
76 RCL 10
77 -
78 RCL 07
79 XEQ "L"
80 ST- 12
81 RCL 10
82 ST+ X
83 ST* 12
84 RCL 09
85 RCL 08
86 RCL 07
87 RCL 10
88 ST+ T
89 +
90 XEQ "L"
91 ST- 12
92 RCL 09
93 RCL 08
94 RCL 07
95 RCL 10
96 ST- T
97 -
98 XEQ "L"
99 ST- 12
100 RCL 09
101 RCL 08
102 RCL 07

103 RCL 10
104 ST+ T
105 -
106 XEQ "L"
107 ST+ 12
108 RCL 09
109 RCL 08
110 RCL 07
111 RCL 10
112 ST- T
113 +
114 XEQ "L"
115 ST+ 12
116 RCL 09
117 RCL 08
118 RCL 10
119 ST+ Z
120 +
121 RCL 07
122 XEQ "L"
123 STO 13
124 RCL 09
125 RCL 08
126 RCL 10
127 ST- Z
128 -
129 RCL 07
130 XEQ "L"
131 ST+ 13
132 RCL 09
133 RCL 08
134 RCL 10
135 ST+ Z
136 -

137 RCL 07
138 XEQ "L"
139 ST- 13
140 RCL 09
141 RCL 08
142 RCL 10
143 ST- Z
144 +
145 RCL 07
146 XEQ "L"
147 RCL 13
148 -
149 RCL 09
150 *
151 ST+ 12
152 RCL 09
153 RCL 10
154 +
155 RCL 08
156 RCL 07
157 XEQ "L"
158 STO 13
159 RCL 09
160 RCL 10
161 -
162 RCL 08
163 RCL 07
164 XEQ "L"
165 ST+ 13
166 RCL 09
167 RCL 08
168 RCL 07
169 XEQ "L"
170 STO 17

171 ST+ X
172 ST- 13
173 RCL 13
174 4
175 *
176 ST/ 12
177 RCL 12
178 RCL 10
179 *
180 2
181 /
182 RCL 09
183 +
184 RCL 10
185 *
186 ST+ 08
187 RCL 12
188 RCL 10
189 ST+ 07
190 *
191 ST+ 09
192 FS? 10
193 GTO 03
194 RCL 17
195 RCL 15
196 X<> 16
197 STO 15
198 *
199 ST+ 14
200 ISG 19
201 CLX
202 RCL 08
203 FC? 00
204 STO IND 19

205 LBL 03
206 DSE 18
207 GTO 03
208 RCL 08
209 RCL 04
210 -
211 FS? 10
212 RTN
213 STO 11
214 RCL 09
215 RCL 08
216 RCL 07
217 XEQ "L"
218 RCL 14
219 +
220 RCL 10
221 *
222 3
223 /
224 STO 07
225 RCL 19
226 .1
227 %
228 20
229 +
230 FS? 00
231 CLX
232 STO 08
233 RCL 09
234 STO 06
235 RCL 05
236 RCL 07
237 CLD
238 END

( 338 bytes / SIZE 021 or 022 + N )

STACK	INPUTS	OUTPUTS
T	/	20.eee or 0
Z	/	y'(b)
Y	m	y'(a)
X	n	I

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:

01 LBL "L"
02 X^2
03 X<>Y
04 X^2
05 +
06 X<>Y
07 X^2
08 +
09 RTN

-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:

x	1	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9	2
y	1	0.9584	0.9266	0.9042	0.8909	0.8867	0.8913	0.9048	0.9273	0.9589	1

-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)

0.9587

0.9270

0.9046

0.8913

0.8868

0.8913

0.9046

0.9270

0.9587

-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

01 LBL "EULAG"
02 STO 06
03 X<>Y
04 STO 05
05 RCL 00
06 2
07 MOD
08 ST+ 00
09 RCL 02
10 RCL 01
11 -
12 RCL 00
13 ST+ X
14 /
15 STO 10
16 SF 10
17 RCL 06
18 XEQ 02
19 STO 14
20 LBL 01
21 VIEW 05
22 RCL 05
23 XEQ 02
24 ENTER
25 ENTER

26 X<> 14
27 -
28 X#0?
29 /
30 RCL 06
31 RCL 05
32 STO 06
33 STO T
34 -
35 *
36 +
37 STO 05
38 -
39 RND
40 X#0?
41 GTO 01
42 2
43 STO 15
44 X^2
45 STO 18
46 20
47 STO 19
48 RCL 05
49 CF 10
50 LBL 02

51 STO 09
52 RCL 01
53 STO 07
54 RCL 00
55 STO 16
56 RCL 03
57 STO 08
58 FS? 10
59 GTO 03
60 RCL 09
61 X<>Y
62 STO 20
63 RCL 07
64 XEQ "L"
65 STO 17
66 LBL 03
67 RCL 09
68 RCL 08
69 RCL 07
70 XEQ "T"
71 RCL 10
72 ST+ 07
73 *
74 STO 12
75 RCL 09

76 +
77 STO 13
78 LASTX
79 RCL 10
80 *
81 STO 11
82 RCL 08
83 +
84 RCL 07
85 XEQ "T"
86 RCL 10
87 ST* 13
88 *
89 ST+ 12
90 ST+ 12
91 RCL 09
92 +
93 ENTER
94 X<> 13
95 ST+ 11
96 ST+ 11
97 RCL 08
98 +
99 RCL 07
100 XEQ "T"

101 RCL 10
102 ST+ 07
103 ST+ X
104 ST* 13
105 *
106 ST+ 12
107 RCL 09
108 +
109 ENTER
110 X<> 13
111 ST+ 11
112 RCL 08
113 +
114 RCL 07
115 XEQ "T"
116 RCL 10
117 ST* 13
118 *
119 RCL 12
120 +
121 3
122 /
123 ST+ 09
124 RCL 11
125 RCL 13

126 +
127 3
128 /
129 ST+ 08
130 FS? 10
131 GTO 03
132 ISG 19
133 CLX
134 RCL 09
135 RCL 08
136 FC? 00
137 STO IND 19
138 RCL 07
139 XEQ "L"
140 RCL 15
141 X<> 18
142 STO 15
143 *
144 ST+ 17
145 LBL 03
146 DSE 16
147 GTO 03
148 RCL 08
149 RCL 04
150 -
151 FS? 10

152 RTN
153 STO 14
154 CLX
155 2
156 /
157 ST- 17
158 RCL 10
159 RCL 17
160 *
161 1.5
162 /
163 STO 07
164 RCL 19
165 .1
166 %
167 20
168 +
169 FS? 00
170 CLX
171 STO 08
172 RCL 09
173 STO 06
174 RCL 05
175 RCL 07
176 CLD
177 END

( 254 bytes / SIZE 021 or 022 + N )

STACK	INPUTS	OUTPUTS
T	/	20.eee or 0
Z	/	y'(b)
Y	m	y'(a)
X	n	I

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:

01 LBL "L"
02 X^2
03 X<>Y
04 X^2
05 +
06 X<>Y
07 X^2
08 +
09 RTN
10 LBL "T"
11 X<>Y
12 RTN

-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:

x	1	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9	2
y	1	0.958715	0.927026	0.904615	0.891257	0.886819	0.891257	0.904615	0.927026	0.958715	1

-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

01 LBL "2ELG"
02 STO 10
03 RDN
04 STO 09
05 RDN
06 STO 08
07 X<>Y
08 STO 07
09 RCL 00
10 2
11 MOD
12 ST+ 00
13 RCL 02
14 RCL 01
15 -
16 RCL 00
17 /
18 STO 15
19 SF 10
20 LBL 01
21 VIEW 07
22 RCL 09
23 RCL 08
24 XEQ 12
25 STO 20
26 X<>Y
27 STO 19
28 RCL 10
29 RCL 07
30 XEQ 12
31 STO 22
32 X<>Y
33 STO 21
34 RCL 09
35 RCL 07
36 XEQ 12
37 STO 18
38 ST- 20
39 ST- 22
40 X<>Y
41 STO 17
42 ST- 19
43 ST- 21
44 RCL 19
45 RCL 22
46 *
47 RCL 20
48 RCL 21
49 *
50 -
51 STO 23
52 X=0?
53 GTO 13
54 RCL 18
55 RCL 21
56 *
57 RCL 17
58 RCL 22
59 *
60 -
61 RCL 23
62 /
63 RCL 08
64 RCL 07
65 STO 08
66 -
67 *
68 ST+ 07
69 ABS
70 RCL 17
71 RCL 20
72 *
73 RCL 18
74 RCL 19
75 *
76 -
77 RCL 23
78 /
79 RCL 10
80 RCL 09
81 STO 10
82 -
83 *
84 ST+ 09
85 ABS
86 +
87 RND
88 X#0?
89 GTO 01
90 LBL 13
91 33
92 STO 32
93 2
94 STO 25
95 X^2

96 STO 24
97 RCL 09
98 RCL 07
99 CF 10
100 GTO 12
101 LBL 07
102 RCL 14
103 RCL 13
104 RCL 12
105 RCL 15
106 ST+ Z
107 +
108 RCL 11
109 XEQ 04
110 STO 30
111 RCL 14
112 RCL 13
113 RCL 12
114 RCL 15
115 ST- Z
116 +
117 RCL 11
118 XEQ 04
119 ST- 30
120 RCL 14
121 RCL 13
122 RCL 12
123 RCL 15
124 ST- Z
125 -
126 RCL 11
127 XEQ 04
128 ST+ 30
129 RCL 14
130 RCL 13
131 RCL 12
132 RCL 15
133 ST+ Z
134 -
135 RCL 11
136 XEQ 04
137 ST- 30
138 RCL 30
139 4
140 /
141 RTN
142 LBL 08
143 RCL 14
144 RCL 13
145 RCL 12
146 RCL 15
147 ST+ T
148 +
149 RCL 11
150 XEQ 05
151 STO 30
152 STO 31
153 RCL 14
154 RCL 13
155 RCL 12
156 RCL 15
157 ST- T
158 -
159 RCL 11
160 XEQ 05
161 ST- 30
162 ST- 31
163 RCL 14
164 RCL 13
165 RCL 12
166 RCL 15
167 ST+ T
168 -
169 RCL 11
170 XEQ 05
171 ST+ 30
172 ST- 31
173 RCL 14
174 RCL 13
175 RCL 12
176 RCL 15
177 ST- T
178 +
179 RCL 11
180 XEQ 05
181 ST- 30
182 ST+ 31
183 RCL 14
184 RCL 15
185 +
186 RCL 13
187 RCL 12
188 RCL 11
189 XEQ 05
190 ST+ X

191 ST- 30
192 RCL 14
193 RCL 15
194 -
195 RCL 13
196 RCL 12
197 RCL 11
198 XEQ 05
199 ST+ X
200 ST+ 30
201 RCL 14
202 RCL 13
203 RCL 12
204 RCL 15
205 +
206 RCL 11
207 XEQ 05
208 ST+ X
209 ST- 31
210 RCL 14
211 RCL 13
212 RCL 12
213 RCL 15
214 -
215 RCL 11
216 XEQ 05
217 ST+ X
218 ST+ 31
219 RCL 31
220 RCL 30
221 2
222 ST/ Z
223 /
224 RTN
225 LBL 09
226 RCL 14
227 RCL 13
228 RCL 12
229 RCL 15
230 ST+ T
231 ST+ Z
232 +
233 RCL 11
234 XEQ 06
235 STO 30
236 RCL 14
237 RCL 13
238 RCL 12
239 RCL 15
240 ST+ T
241 ST- Z
242 +
243 RCL 11
244 XEQ 06
245 ST- 30
246 RCL 14
247 RCL 13
248 RCL 12
249 RCL 15
250 ST+ T
251 ST- Z
252 -
253 RCL 11
254 XEQ 06
255 ST+ 30
256 RCL 14
257 RCL 13
258 RCL 12
259 RCL 15
260 ST+ T
261 ST+ Z
262 -
263 RCL 11
264 XEQ 06
265 ST- 30
266 RCL 14
267 RCL 13
268 RCL 12
269 RCL 15
270 ST- T
271 ST+ Z
272 +
273 RCL 11
274 XEQ 06
275 ST- 30
276 RCL 14
277 RCL 13
278 RCL 12
279 RCL 15
280 ST- T
281 ST- Z
282 +
283 RCL 11
284 XEQ 06
285 ST+ 30

286 RCL 14
287 RCL 13
288 RCL 12
289 RCL 15
290 ST- T
291 ST- Z
292 -
293 RCL 11
294 XEQ 06
295 ST- 30
296 RCL 14
297 RCL 13
298 RCL 12
299 RCL 15
300 ST- T
301 ST+ Z
302 -
303 RCL 11
304 XEQ 06
305 RCL 30
306 +
307 RTN
308 LBL 04
309 FS? 02
310 X<>Y
311 FC? 03
312 GTO "L"
313 R^
314 X<> T
315 RDN
316 GTO "L"
317 LBL 05
318 FS? 02
319 X<>Y
320 FC? 03
321 GTO "L"
322 RDN
323 X<>Y
324 R^
325 GTO "L"
326 LBL 06
327 FS? 02
328 X<>Y
329 FC? 03
330 GTO "L"
331 X<> Z
332 X<>Y
333 GTO "L"
334 LBL 12
335 STO 14
336 X<>Y
337 STO 16
338 RCL 05
339 STO 13
340 RCL 00
341 STO 26
342 RCL 03
343 STO 12
344 RCL 01
345 STO 11
346 FS? 10
347 GTO 03
348 RCL 14
349 RCL 13
350 RCL 12
351 RCL 11
352 XEQ "L"
353 CHS
354 STO 27
355 LBL 03
356 RCL 14
357 RCL 13
358 RCL 15
359 +
360 RCL 12
361 RCL 11
362 XEQ "L"
363 STO 29
364 RCL 14
365 RCL 13
366 RCL 12
367 RCL 11
368 XEQ "L"
369 STO 28
370 ST+ X
371 ST- 29
372 RCL 14
373 RCL 13
374 RCL 15
375 -
376 RCL 12
377 RCL 11
378 XEQ "L"
379 RCL 29
380 +

381 RCL 14
382 *
383 STO 29
384 RCL 14
385 RCL 13
386 RCL 12
387 RCL 15
388 +
389 RCL 11
390 XEQ "L"
391 STO 30
392 RCL 14
393 RCL 13
394 RCL 12
395 RCL 15
396 -
397 RCL 11
398 XEQ "L"
399 ST- 30
400 RCL 30
401 RCL 15
402 2
403 /
404 *
405 ST- 29
406 CF 02
407 CF 03
408 XEQ 07
409 RCL 13
410 *
411 ST+ 29
412 RCL 11
413 X<> 12
414 STO 11
415 SF 02
416 XEQ 07
417 ST+ 29
418 RCL 11
419 X<> 12
420 STO 11
421 RCL 13
422 X<> 14
423 STO 13
424 CF 02
425 SF 03
426 XEQ 07
427 RCL 13
428 *
429 ST- 29
430 RCL 15
431 ST* 29
432 RCL 13
433 X<> 14
434 STO 13
435 CF 03
436 XEQ 08
437 RCL 13
438 ST* Z
439 X^2
440 *
441 ST- 29
442 CLX
443 RCL 16
444 *
445 ST+ X
446 ST- 29
447 RCL 11
448 X<> 12
449 STO 11
450 SF 02
451 XEQ 08
452 ST- 29
453 CLX
454 RCL 16
455 ST+ X
456 *
457 ST- 29
458 RCL 11
459 X<> 12
460 STO 11
461 RCL 12
462 X<> 13
463 STO 12
464 CF 02
465 SF 03
466 XEQ 08
467 RCL 14
468 ST* Z
469 X^2
470 *
471 ST- 29
472 CLX
473 RCL 16
474 ST+ X
475 *

476 ST- 29
477 RCL 12
478 X<> 13
479 STO 12
480 CF 03
481 XEQ 09
482 RCL 13
483 *
484 RCL 14
485 *
486 ST- 29
487 RCL 11
488 X<> 12
489 STO 11
490 SF 02
491 XEQ 09
492 RCL 14
493 *
494 ST- 29
495 RCL 11
496 X<> 13
497 STO 11
498 CF 02
499 SF 03
500 XEQ 09
501 RCL 11
502 *
503 ST- 29
504 RCL 11
505 X<> 13
506 X<> 12
507 STO 11
508 RCL 14
509 RCL 15
510 ST+ X
511 +
512 RCL 13
513 RCL 12
514 RCL 11
515 CF 03
516 XEQ "L"
517 2
518 /
519 STO 30
520 RCL 14
521 RCL 15
522 +
523 RCL 13
524 RCL 12
525 RCL 11
526 XEQ "L"
527 ST- 30
528 STO 31
529 RCL 14
530 RCL 15
531 -
532 RCL 13
533 RCL 12
534 RCL 11
535 XEQ "L"
536 ST+ 30
537 ST+ 31
538 RCL 14
539 RCL 15
540 ST+ X
541 -
542 RCL 13
543 RCL 12
544 RCL 11
545 XEQ "L"
546 2
547 /
548 ST- 30
549 RCL 30
550 RCL 16
551 X^2
552 *
553 ST- 29
554 RCL 29
555 RCL 31
556 RCL 28
557 ST+ X
558 -
559 RCL 15
560 *
561 /
562 STO 29
563 RCL 15
564 *
565 4
566 /
567 RCL 16
568 +
569 RCL 15
570 *

571 3
572 /
573 RCL 14
574 +
575 RCL 15
576 *
577 2
578 /
579 RCL 13
580 +
581 RCL 15
582 *
583 ST+ 12
584 RCL 29
585 RCL 15
586 *
587 3
588 /
589 RCL 16
590 +
591 RCL 15
592 *
593 2
594 /
595 RCL 14
596 +
597 RCL 15
598 *
599 ST+ 13
600 RCL 29
601 RCL 15
602 *
603 2
604 /
605 RCL 16
606 +
607 RCL 15
608 *
609 ST+ 14
610 RCL 29
611 RCL 15
612 ST+ 11
613 *
614 ST+ 16
615 FS? 10
616 GTO 03
617 RCL 28
618 RCL 24
619 X<> 25
620 STO 24
621 *
622 ST+ 27
623 ISG 32
624 CLX
625 RCL 12
626 FC? 00
627 STO IND 32
628 LBL 03
629 DSE 26
630 GTO 03
631 RCL 13
632 RCL 06
633 -
634 RCL 12
635 RCL 04
636 -
637 FS? 10
638 RTN
639 RCL 14
640 RCL 13
641 RCL 12
642 RCL 11
643 XEQ "L"
644 RCL 27
645 +
646 RCL 15
647 *
648 3
649 /
650 STO 11
651 RCL 16
652 STO 10
653 RCL 32
654 .1
655 %
656 33
657 +
658 FS? 00
659 CLX
660 STO 12
661 RCL 14
662 STO 08
663 RCL 07
664 RCL 11
665 CLD
666 END

( 985 bytes SIZE 033 or 34+N )

STACK	INPUTS	OUTPUTS
T	m	33.eee or 0
Z	n	y''(b)
Y	m'	y''(a)
X	n'	I

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:

01 LBL "L"
02 X<>Y
03 X^2
04 +
05 X<>Y
06 X^2
07 ST+ X
08 +
09 X<>Y
10 X^2
11 +
12 RTN

-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:

x	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
y	0	0.1004	0.2037	0.3131	0.4318	0.5631	0.7106	0.8781	1.0696	1.2897	1.5431

-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

0.1005

0.2040

0.3136

0.4324

0.5638

0.7113

0.8786

1.0699

1.2898

1.5431

-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

01 LBL "2ELG"
02 STO 10
03 RDN
04 STO 09
05 RDN
06 STO 08
07 X<>Y
08 STO 07
09 RCL 00
10 2
11 MOD
12 ST+ 00
13 RCL 02
14 RCL 01
15 -
16 RCL 00
17 ST+ X
18 /
19 STO 16
20 SF 10
21 LBL 01
22 VIEW 07
23 RCL 09
24 RCL 08
25 XEQ 12
26 STO 33
27 X<>Y
28 STO 32
29 RCL 10
30 RCL 07
31 XEQ 12
32 STO 35
33 X<>Y
34 STO 34
35 RCL 09
36 RCL 07
37 XEQ 12
38 STO 31
39 ST- 33
40 ST- 35
41 X<>Y
42 STO 30
43 ST- 32

44 ST- 34
45 RCL 32
46 RCL 35
47 *
48 RCL 33
49 RCL 34
50 *
51 -
52 STO 36
53 X=0?
54 GTO 13
55 RCL 31
56 RCL 34
57 *
58 RCL 30
59 RCL 35
60 *
61 -
62 RCL 36
63 /
64 RCL 08
65 RCL 07
66 STO 08
67 -
68 *
69 ST+ 07
70 ABS
71 RCL 30
72 RCL 33
73 *
74 RCL 31
75 RCL 32
76 *
77 -
78 RCL 36
79 /
80 RCL 10
81 RCL 09
82 STO 10
83 -
84 *
85 ST+ 09
86 ABS

87 +
88 RND
89 X#0?
90 GTO 01
91 LBL 13
92 37
93 STO 29
94 2
95 STO 26
96 X^2
97 STO 27
98 RCL 09
99 RCL 07
100 CF 10
101 LBL 12
102 STO 14
103 X<>Y
104 STO 15
105 RCL 05
106 STO 13
107 RCL 00
108 STO 25
109 RCL 01
110 STO 11
111 RCL 03
112 STO 12
113 FS? 10
114 GTO 03
115 STO 37
116 XEQ "L"
117 STO 28
118 LBL 03
119 XEQ "T"
120 RCL 16
121 ST+ 11
122 *
123 STO 24
124 RCL 15
125 STO 20
126 +
127 STO 15
128 LASTX
129 RCL 16

130 *
131 STO 23
132 RCL 14
133 STO 19
134 +
135 STO 14
136 LASTX
137 RCL 16
138 *
139 STO 22
140 RCL 13
141 STO 18
142 +
143 STO 13
144 LASTX
145 RCL 16
146 *
147 STO 21
148 RCL 12
149 STO 17
150 +
151 STO 12
152 XEQ "T"
153 RCL 16
154 *
155 ST+ 24
156 ST+ 24
157 RCL 20
158 +
159 X<> 15
160 RCL 16
161 *
162 ST+ 23
163 ST+ 23
164 RCL 19
165 +
166 X<> 14
167 RCL 16
168 *
169 ST+ 22
170 ST+ 22
171 RCL 18
172 +

173 X<> 13
174 RCL 16
175 *
176 ST+ 21
177 ST+ 21
178 RCL 17
179 +
180 STO 12
181 XEQ "T"
182 RCL 16
183 ST+ 11
184 ST+ X
185 *
186 ST+ 24
187 RCL 20
188 +
189 X<> 15
190 RCL 16
191 ST+ X
192 *
193 ST+ 23
194 RCL 19
195 +
196 X<> 14
197 RCL 16
198 ST+ X
199 *
200 ST+ 22
201 RCL 18
202 +
203 X<> 13
204 RCL 16
205 ST+ X
206 *
207 ST+ 21
208 RCL 17
209 +
210 STO 12
211 XEQ "T"
212 RCL 16
213 *
214 ST+ 24
215 RCL 15

216 LASTX
217 *
218 ST+ 23
219 RCL 14
220 LASTX
221 *
222 ST+ 22
223 RCL 13
224 LASTX
225 *
226 ST+ 21
227 RCL 21
228 3
229 /
230 RCL 17
231 +
232 STO 12
233 RCL 22
234 3
235 /
236 RCL 18
237 +
238 STO 13
239 RCL 23
240 3
241 /
242 RCL 19
243 +
244 STO 14
245 RCL 24
246 3
247 /
248 RCL 20
249 +
250 STO 15
251 FS? 10
252 GTO 03
253 ISG 29
254 CLX
255 RCL 12
256 FC? 00
257 STO IND 29
258 XEQ "L"

259 RCL 26
260 X<> 27
261 STO 26
262 *
263 ST+ 28
264 LBL 03
265 DSE 25
266 GTO 03
267 RCL 13
268 RCL 06
269 -
270 RCL 12
271 RCL 04
272 -
273 FS? 10
274 RTN
275 R^
276 2
277 /
278 ST- 28
279 RCL 28
280 RCL 16
281 *
282 1.5
283 /
284 STO 11
285 RCL 15
286 STO 10
287 RCL 14
288 STO 08
289 RCL 29
290 .1
291 %
292 37
293 +
294 FS? 00
295 CLX
296 STO 12
297 RCL 09
298 RCL 07
299 RCL 11
300 CLD
301 END

( 460 bytes SIZE 037 or 38+N )

STACK	INPUTS	OUTPUTS
T	m	37.eee or 0
Z	n	y'''(a)
Y	m'	y''(a)
X	n'	I

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:

01 LBL "L"
02 RCL 11
03 RCL 12
04 X^2
05 +
06 RCL 13
07 X^2
08 ST+ X
09 +
10 RCL 14
11 X^2
12 +
13 RTN
14 LBL "T"
15 RCL 14
16 ST+ X
17 RCL 12
18 -
19 RTN

-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:

x	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
y	0	0.100500	0.204013	0.313601	0.432429	0.563813	0.711280	0.878619	1.069949	1.289778	1.543081

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

0.100500

0.204013

0.313602

0.432429

0.563813

0.711279

0.878618

1.069948

1.289778

1.543081

-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

01 LBL "ELG2"
02 STO 10
03 RDN
04 STO 09
05 RDN
06 STO 08
07 X<>Y
08 STO 07
09 RCL 00
10 2
11 MOD
12 ST+ 00
13 RCL 02
14 RCL 01
15 -
16 RCL 00
17 /
18 STO 16
19 ST+ X
20 STO 17
21 11
22 STO 20
23 12
24 STO 21
25 13
26 STO 22
27 14
28 STO 23
29 15
30 STO 24
31 SF 10
32 LBL 01
33 VIEW 07
34 RCL 09
35 RCL 08
36 XEQ 12
37 STO 33
38 X<>Y
39 STO 32
40 RCL 10
41 RCL 07
42 XEQ 12
43 STO 35
44 X<>Y
45 STO 34
46 RCL 09
47 RCL 07
48 XEQ 12
49 STO 31
50 ST- 33
51 ST- 35

52 X<>Y
53 STO 30
54 ST- 32
55 ST- 34
56 RCL 32
57 RCL 35
58 *
59 RCL 33
60 RCL 34
61 *
62 -
63 STO 36
64 X=0?
65 GTO 13
66 RCL 31
67 RCL 34
68 *
69 RCL 30
70 RCL 35
71 *
72 -
73 RCL 36
74 /
75 RCL 08
76 RCL 07
77 STO 08
78 -
79 *
80 ST+ 07
81 ABS
82 RCL 30
83 RCL 33
84 *
85 RCL 31
86 RCL 32
87 *
88 -
89 RCL 36
90 /
91 RCL 10
92 RCL 09
93 STO 10
94 -
95 *
96 ST+ 09
97 ABS
98 +
99 RND
100 X#0?
101 GTO 01
102 LBL 13

103 43
104 STO 38
105 2
106 STO 39
107 X^2
108 STO 40
109 RCL 09
110 RCL 07
111 CF 10
112 GTO 12
113 LBL 07
114 RCL 16
115 ST+ IND 25
116 ST+ IND 26
117 XEQ "L"
118 STO 27
119 RCL 17
120 ST- IND 26
121 XEQ "L"
122 ST- 27
123 RCL 17
124 ST- IND 25
125 XEQ "L"
126 ST+ 27
127 RCL 17
128 ST+ IND 26
129 XEQ "L"
130 ST- 27
131 RCL 16
132 ST+ IND 25
133 ST- IND 26
134 RCL 27
135 RTN
136 LBL 12
137 STO 14
138 X<>Y
139 STO 15
140 RCL 05
141 STO 13
142 RCL 00
143 STO 41
144 RCL 01
145 STO 11
146 RCL 03
147 STO 12
148 FS? 10
149 GTO 03
150 STO 43
151 44
152 RCL 00
153 +

154 RCL 13
155 FC? 00
156 STO IND Y
157 XEQ "L"
158 CHS
159 STO 42
160 LBL 03
161 RCL 16
162 ST+ 12
163 XEQ "L"
164 STO 18
165 RCL 17
166 ST- 12
167 XEQ "L"
168 ST- 18
169 RCL 16
170 ST+ 12
171 ST+ 13
172 XEQ "L"
173 STO 19
174 RCL 17
175 ST- 13
176 XEQ "L"
177 ST- 19
178 RCL 16
179 ST+ 13
180 RCL 17
181 ST* 18
182 ST* 19
183 RCL 20
184 STO 25
185 RCL 23
186 STO 26
187 XEQ 07
188 ST- 18
189 RCL 24
190 STO 26
191 XEQ 07
192 ST- 19
193 RCL 21
194 STO 25
195 XEQ 07
196 RCL 14
197 *
198 ST- 19
199 RCL 23
200 STO 26
201 XEQ 07
202 RCL 14
203 *
204 ST- 18

205 RCL 22
206 STO 25
207 XEQ 07
208 RCL 15
209 *
210 ST- 18
211 RCL 24
212 STO 26
213 XEQ 07
214 RCL 15
215 *
216 ST- 19
217 RCL 23
218 STO 25
219 XEQ 07
220 STO 28
221 RCL 16
222 ST+ 14
223 XEQ "L"
224 STO 29
225 RCL 17
226 ST- 14
227 XEQ "L"
228 ST+ 29
229 RCL 16
230 ST+ 14
231 ST+ 15
232 XEQ "L"
233 STO 30
234 RCL 17
235 ST- 15
236 XEQ "L"
237 ST+ 30
238 RCL 16
239 ST+ 15
240 XEQ "L"
241 STO 37
242 ST+ X
243 ST- 29
244 ST- 30
245 4
246 ST* 29
247 ST* 30
248 RCL 18
249 RCL 30
250 *
251 RCL 19
252 RCL 28
253 *
254 -
255 RCL 29

256 RCL 30
257 *
258 RCL 28
259 X^2
260 -
261 STO 26
262 /
263 STO 25
264 RCL 19
265 RCL 29
266 *
267 RCL 18
268 RCL 28
269 *
270 -
271 RCL 26
272 /
273 STO 26
274 RCL 16
275 *
276 2
277 /
278 RCL 15
279 +
280 RCL 16
281 *
282 ST+ 13
283 RCL 25
284 RCL 16
285 *
286 2
287 /
288 RCL 14
289 +
290 RCL 16
291 *
292 ST+ 12
293 RCL 26
294 LASTX
295 *
296 ST+ 15
297 RCL 25
298 LASTX
299 ST+ 11
300 *
301 ST+ 14
302 FS? 10
303 GTO 03
304 ISG 38
305 CLX
306 RCL 12

307 FC? 00
308 STO IND 38
309 RCL 38
310 RCL 00
311 +
312 1
313 +
314 RCL 13
315 FC? 00
316 STO IND Y
317 RCL 37
318 RCL 40
319 X<> 39
320 STO 40
321 *
322 ST+ 42
323 LBL 03
324 DSE 41
325 GTO 03
326 RCL 13
327 RCL 06
328 -
329 RCL 12
330 RCL 04
331 -
332 FS? 10
333 RTN
334 XEQ "L"
335 RCL 42
336 +
337 RCL 16
338 *
339 3
340 /
341 STO 11
342 RCL 15
343 STO 10
344 RCL 14
345 STO 08
346 RCL 38
347 .1
348 %
349 43
350 +
351 FS? 00
352 CLX
353 STO 12
354 RCL 09
355 RCL 07
356 RCL 11
357 CLD
358 END

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

STACK	INPUTS	OUTPUTS
T	m	43.eee or 0
Z	n	z'(a)
Y	m'	y'(a)
X	n'	I

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:

01 LBL "L"
02 RCL 12
03 X^2
04 RCL 11
05 RCL 13
06 *
07 ST* Y
08 X^2
09 -
10 RCL 14
11 X^2
12 +
13 RCL 15
14 X^2
15 2
16 /
17 +
18 RTN

-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

x	1	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9	2
y	1	1.0786	1.1575	1.2374	1.3192	1.4047	1.4960	1.5963	1.7098	1.8420	2.0000
z	0	0.0436	0.0981	0.1651	0.2459	0.3412	0.4510	0.5743	0.7090	0.8521	1.0000

-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

01 LBL "ELG2"
02 STO 10
03 RDN
04 STO 09
05 RDN
06 STO 08
07 X<>Y
08 STO 07
09 RCL 00
10 2
11 MOD
12 ST+ 00
13 RCL 02
14 RCL 01
15 -
16 RCL 00
17 ST+ X
18 /
19 STO 16
20 SF 10
21 LBL 01
22 VIEW 07
23 RCL 09
24 RCL 08
25 XEQ 12
26 STO 33
27 X<>Y
28 STO 32
29 RCL 10
30 RCL 07
31 XEQ 12
32 STO 35
33 X<>Y
34 STO 34
35 RCL 09
36 RCL 07
37 XEQ 12
38 STO 31
39 ST- 33
40 ST- 35
41 X<>Y
42 STO 30
43 ST- 32
44 ST- 34
45 RCL 32

46 RCL 35
47 *
48 RCL 33
49 RCL 34
50 *
51 -
52 STO 36
53 X=0?
54 GTO 13
55 RCL 31
56 RCL 34
57 *
58 RCL 30
59 RCL 35
60 *
61 -
62 RCL 36
63 /
64 RCL 08
65 RCL 07
66 STO 08
67 -
68 *
69 ST+ 07
70 ABS
71 RCL 30
72 RCL 33
73 *
74 RCL 31
75 RCL 32
76 *
77 -
78 RCL 36
79 /
80 RCL 10
81 RCL 09
82 STO 10
83 -
84 *
85 ST+ 09
86 ABS
87 +
88 RND
89 X#0?
90 GTO 01

91 LBL 13
92 37
93 STO 29
94 2
95 STO 26
96 X^2
97 STO 27
98 RCL 09
99 RCL 07
100 CF 10
101 LBL 12
102 STO 14
103 X<>Y
104 STO 15
105 RCL 05
106 STO 13
107 RCL 00
108 STO 25
109 RCL 01
110 STO 11
111 RCL 03
112 STO 12
113 FS? 10
114 GTO 03
115 STO 37
116 38
117 RCL 00
118 +
119 RCL 13
120 FC? 00
121 STO IND Y
122 XEQ "L"
123 STO 28
124 LBL 03
125 XEQ "T"
126 RCL 16
127 ST+ 11
128 ST* Z
129 *
130 STO 24
131 RCL 15
132 STO 20
133 +
134 STO 15
135 LASTX

136 RCL 16
137 *
138 STO 22
139 RCL 13
140 STO 18
141 +
142 STO 13
143 X<> L
144 RCL 16
145 *
146 STO 21
147 R^
148 STO 23
149 RCL 14
150 STO 19
151 +
152 STO 14
153 LASTX
154 RCL 16
155 *
156 STO 21
157 RCL 12
158 STO 17
159 +
160 STO 12
161 XEQ "T"
162 RCL 16
163 ST* Z
164 *
165 ST+ 24
166 ST+ 24
167 RCL 20
168 +
169 X<> 15
170 RCL 16
171 *
172 ST+ 22
173 ST+ 22
174 RCL 18
175 +
176 STO 13
177 X<>Y
178 ST+ 23
179 ST+ 23
180 RCL 19

181 +
182 X<> 14
183 RCL 16
184 *
185 ST+ 21
186 ST+ 21
187 RCL 17
188 +
189 STO 12
190 XEQ "T"
191 RCL 16
192 ST+ 11
193 ST+ X
194 ST* Z
195 *
196 ST+ 24
197 RCL 20
198 +
199 X<> 15
200 RCL 16
201 ST+ X
202 *
203 ST+ 22
204 RCL 18
205 +
206 STO 13
207 X<>Y
208 ST+ 23
209 RCL 19
210 +
211 X<> 14
212 RCL 16
213 ST+ X
214 *
215 ST+ 21
216 RCL 17
217 +
218 STO 12
219 XEQ "T"
220 RCL 16
221 ST* Z
222 *
223 ST+ 24
224 RCL 15
225 LASTX

226 *
227 ST+ 22
228 R^
229 ST+ 23
230 RCL 14
231 LASTX
232 *
233 ST+ 21
234 3
235 ST/ 21
236 ST/ 22
237 ST/ 23
238 ST/ 24
239 RCL 17
240 RCL 21
241 +
242 STO 12
243 RCL 18
244 RCL 22
245 +
246 STO 13
247 RCL 19
248 RCL 23
249 +
250 STO 14
251 RCL 20
252 RCL 24
253 +
254 STO 15
255 FS? 10
256 GTO 03
257 ISG 29
258 CLX
259 RCL 12
260 FC? 00
261 STO IND 29
262 RCL 29
263 RCL 00
264 +
265 1
266 +
267 RCL 13
268 FC? 00
269 STO IND Y

270 XEQ "L"
271 RCL 26
272 X<> 27
273 STO 26
274 *
275 ST+ 28
276 LBL 03
277 DSE 25
278 GTO 03
279 RCL 13
280 RCL 06
281 -
282 RCL 12
283 RCL 04
284 -
285 FS? 10
286 RTN
287 R^
288 2
289 /
290 ST- 28
291 RCL 28
292 RCL 16
293 *
294 1.5
295 /
296 STO 11
297 RCL 15
298 STO 10
299 RCL 14
300 STO 08
301 RCL 29
302 .1
303 %
304 37
305 +
306 FS? 00
307 CLX
308 STO 12
309 RCL 09
310 RCL 07
311 RCL 11
312 CLD
313 END

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

STACK	INPUTS	OUTPUTS
T	m	37.eee or 0
Z	n	z'(a)
Y	m'	y'(a)
X	n'	I

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:

01 LBL "L"
02 RCL 12
03 X^2
04 RCL 11
05 RCL 13
06 *
07 ST* Y
08 X^2
09 -
10 RCL 14
11 X^2
12 +
13 RCL 15
14 X^2
15 2
16 /
17 +
18 RTN
19 LBL "T"
20 RCL 12
21 ENTER^
22 X^2
23 RCL 13
24 ST* Z
25 ST+ X
26 RCL 11
27 ST* T
28 ST*Z
29 X^2
30 *
31 -
32 RTN

-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

x	1	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9	2
y	1	1.0721	1.1448	1.2191	1.2964	1.3788	1.4689	1.5706	1.6887	1.8292	2.0000
z	0	0.0506	0.1126	0.1871	0.2745	0.3745	0.4862	0.6074	0.7354	0.8672	1.0000

-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

hp41programs

Euler-Lagrange Equation for the HP-41