(48G 49G 50G) Runge-Kutta Diff Eqn solver

[rcthacker]

I would like to use the Runge-Kutta-Fehlberg function/solver for second order nonlinear DE's. It works fine for first order NL equations, and for second order linear equations. However, I can't figure out how to use the matrix (function) - vector (solution) multiplication to include nonlinear dependent variables. For example, y' = v, and y'' = g - k X v^2. (the skydiver problem from calculus). BTW, the 50G User Guide discusses first order linear equations, but solves a first order nonlinear equation (SQRT(V)).

Any recommendations for further reading on the use of this powerful function?

[Werner]

You can rewrite that as

[ y v ]' = [ v g-k.X.v^2 ]

So the ODE program must accept [y v] and return [v g-k.X.v^2]
Something like this, with g,k and X global vars:

Code:

\<<
  2. GET @ v
  g OVER SQ
  k * X * - \->NUM
  2. \->ARRY
\>>

Cheers, Werner

[rcthacker]

Thank you, Werner, for your comment.

Could you provide some clarifications? What do \ and @ in your Code mean or do? I don't see them in any programs in the Advanced User Reference. If I key your suggested procedure (between the guillemets) nothing is saved from the @ on. Is [y v] the input to the procedure? ( 2 GET . . .) And, ARRY the return? How is the procedure called? It can't be substituted in the solution vector [y v].

I made a mistake in my question: should be y' = G - K * v^2 (* looked too small for legibility, so I used X). G = 9.8 and K = 0.005. So, y' = 9.8 - 0.005 * v^2.

Thanks for any help.

[BINUBALL]

@ is program comment. Although they will disappear after you enter program. On PC, they are not disappear.
Some command including \ is escaping. 50g program can be made in computer. If you have made some programs on 50g. It needs some character that can't write on PC such as ->, <<, >>. So you write \<< on PC and send to calculator, using ascii to binary converter. And \<< converted into <<.

p.s. so \->ARRY is ->ARRY on 50g

[Werner]

Having never used the differential equation solver in the 49/50, I had to look up how to use it.
Given Y' = F(X,Y), with X the independent and Y the dependent variable, it expects an equation or a program to return the value of Y', given X and Y.
In your case, you have a second-order Diff.Eq
y' = v
v' = g - k*v^2
that we'll rewrite as a first-order one:
U' = F(U,t)
with U = [ y v ]
and F(U,t) = [ v g-k*v^2 ]
Instead of an equation for F, the calculator also accepts a program, so all we need to do
is write one that will return [ v g-k*v^2 ] from U = [ y v ]

Code:

\<<
  U 2. GET
  g OVER SQ k * - \->NUM
  2. \->ARRY
\>>

Seeing that the program is very small, just type it in your 50g, without the backslashes ;-)
Store the program in the reserved variable EQ, fire up the num. solver, diff eq.
The program is shown as 'F', enter indep. var' t' and soln 'U' and all the other parameters, and SOLVE. Should work.

Update: two more remarks
- the equation should be v' = -g - k*v^2 (you're falling down, no?)
- the inital value for U must be given as a 1x2 vector [ y0 v0 ] eg [ 1000. 0. ] so falling from 1000m high, initial velocity 0. m/s
Solving for Final will then give the final value of U [ vf vf' ]

Cheers, Werner

[rcthacker]

Perfect!

Thank you, Werner.

"Store the program in the reserved variable EQ, fire up the num. solver, diff eq.
The program is shown as 'F'," This is the important part that should be included in the 50G documents. It allows the solution of 2nd order nonlinear ODE's. A very important class.

Thanks, also, to Eric Rechlin who sent me to this forum.

rcthacker

[Werner]

You can also use EDIT in the DiffEq form itself to create the program.
Werner

[rcthacker]

"Instead of an equation for F, the calculator also accepts a program, so all we need to do is write one that will return [ v g-k*v^2 ] from U = [ y v ]" This isn't in any documentation I've seen. Where did you find it?

This works fine from num.slv.

What would you do if invoking from the program function RKF? I tried it, but it wouldn't accept the program inlieu of a function.

[Werner]

I tried it ;-)
It helps if you know that equations are really stored as programs.

As for RKF: it works if you store your program in a variable. Since it is already in EQ, that's what I used. So, for the above problem:

0. 't' STO ( t0 )
[ 1000. 0. ] 'U' STO ( U0 )
{ t U EQ }
0.001 (accepting default stepsize )
3. (or anything, tn)
RKF

Cheers, Werner

[rcthacker]

<< U 2 GET 9.8 OVER SQ .005 * - 2 -> ARRY >> 'EQ STO 0 'T STO [0 0] 'U STO

In NUM.SLV F: EQ program Tfinal: 10 U solve: [306..., 43....], expected solution

As a program called function:

With 0 'T STO [0 0] 'U STO

{ T U EQ}
.0001
10
RKF -> Bad Argument Type

What now?

[Werner]

You're in exact mode, and so [0 0] is NOT the same as [0. 0.].
The former is a symbolic matrix, the latter a numeric, and that is what we need here (and why you're getting the 'Bad Argument Type')
Provide values with decimal points throughout and it will work.
(or, switch on approx. mode: MODE/CAS and select Approx).

Cheers, Werner

[rcthacker]

Thank you again, Werner. More thoughts on your responses:

"Having never used the differential equation solver in the 49/50, I had to look up how to use it.
Given Y' = F(X,Y), with X the independent and Y the dependent variable, it expects an equation or a program to return the value of Y',

"Instead of an equation for F, the calculator also accepts a program,

"It helps if you know that equations are really stored as programs.

Equations are stored as programs?

How or where did you learn that a procedure could be substituted for the RHS of the DE? Was this characteristic also in the 48G and 49G firmware? Interesting that the issue was never raised back then. Or, I suppose it would have been addressed in succeeding User Guides/References The work-around or its need does not seem consistent with a straight-forward Runge-Kutta formulation. Somewhat troubling.

[Werner]

(06-11-2021 11:23 PM)rcthacker Wrote:  Equations are stored as programs?

Yes. Equations are 'compiled' into programs, and turned into equations again for stack display. On my 49G, doing
'(A+B)*(C+D)' \->S2
results in

Code:

SYMBOL
 ID A
 ID B
 x*
 ID C
 ID D
 x*
 x+
;
@

which is the way the equation is stored, with SysRPL commands. (anything starting with a small 'x' is a user command, ID is a global name, SYMBOL defines an equation and ; ends it).

Quote:How or where did you learn that a procedure could be substituted for the RHS of the DE?

As I said, I just tried it out. The ODE solver just needs an 'object' to return f(x,y), it does not manipulate the equation. All equations have program equivalents, but not all programs can be turned into equations, so I would have been very surprised if the input had been limited to equations only.

Quote:Was this characteristic also in the 48G and 49G firmware?

I only have a 48 and a 49. They both have the ODE solver, and I tried out things on the 49.

Quote:The work-around or its need does not seem consistent with a straight-forward Runge-Kutta formulation. Somewhat troubling.

I don't understand what you mean here.

Cheers, Werner

[rcthacker]

Again, thank you, Werner, for your enlightening comments and discussion.

I went through my copy of Bill Wickes' HP-28 Insights looking for Expressions, Equations, Programs, and their equivalence. Probably the best description is on page 90, 6.3 Symbolic Math (First Edition, First Printing):

"As we will show in section 9.1, algebraic objects are procedures that are internally the same as programs. So just creating any algebraic object is equivalent to writing a program. Its "inputs" are the values stored in the variables named within the algebraic object; its "output" is the symbolic or numeric result that is returned to the stack. . . "

I couldn't find a copy of his HP-48 Insights, but I imagine it would say the same. In section 9.1 he does say, in effect, that algebraics are more restrictive.

"All equations have program equivalents, but not all programs can be turned into equations, so I would have been very surprised if the input had been limited to equations only." I doubt the firmware writers designed the input arguments to allow for this. Sort of like an unsupported feature. But, maybe it is a central if undescribed characteristic of the RPL design.

Since I was unable to access the necessary variable in my original attempt to solve the 2nd order DE, it was very valuable that you showed how an RPN procedure could provide the required "inputs" and return the required "outputs". All the examples in the 48G and 50G Guides and AURs (nothing for the 49G) show only algebraics or a matrix-vector product - no RPN procedures. The problem could only be solved using RPN commands not allowed in an algebraic object. In a straight-forward Runge-Kutta code all of the function inputs and outputs would be accessible variables. It's a mystery to me that the firmware writers decided to represent 2nd order DEs with a matrix of coefficients and a vector of dependent variables: neither representable in an algebraic expression/object in the RPL OS/User language. Maybe an academic feature.

[Werner]

(06-14-2021 10:24 PM)rcthacker Wrote:  "All equations have program equivalents, but not all programs can be turned into equations, so I would have been very surprised if the input had been limited to equations only." I doubt the firmware writers designed the input arguments to allow for this.

If all they wanted to allow were equations, it would have been very easy to restrict the inputs to equations only, but they didn't.

Quote:In a straight-forward Runge-Kutta code all of the function inputs and outputs would be accessible variables. It's a mystery to me that the firmware writers decided to represent 2nd order DEs with a matrix of coefficients and a vector of dependent variables: neither representable in an algebraic expression/object in the RPL OS/User language. Maybe an academic feature.

The firmware writers didn't decide that; it's a common way of representing higher-order ODEs as first-order ones.

Cheers, Werner

[Giuseppe Donnini]

(06-14-2021 10:24 PM)rcthacker Wrote:  I doubt the firmware writers designed the input arguments to allow for this.

Before commenting on the intentions of the firmware designers, I would first listen to what the very head of the firmware design team, Dr. William Wickes, has to say about it: on page 415 of his HP 48 Insights, Part II: Problem Solving Ressources (of which a free copy is available here), we may read the following:

A nice property of HP 48 Solve is that the current “equation” does not literally have to be an equation. You can use programs, expressions, and equations almost interchangeably for solving purposes. In effect, HP Solve always solves f(x) = 0. In the case of an HP 48 algebraic expression object, f(x) is the expression represented by the object, where x is the name of the unknown variable. For an equation object representing g(x) = h(x), HP Solve solves f(x) = g(x)-h(x) = 0. For a program, f(x) is the expression that is equivalent to that (RPN) program. Programs used for this purpose must be equivalent to an algebraic expression, taking no arguments from the stack and returning one real number or unit object.

And an equivalent, albeit shorter, explanation is already given on page 97 of his HP-28 Insights, which you claim you have read.

So, yes, this is indeed a central characteristic of the RPL design.

[rcthacker]

HP Solve doesn't solve DEs numerically. Runge-Kutta integrates along the unknown function trajectory. HP Solve finds roots of expressions. Not the same. I was interested in Vfinal and Sfinal when
V' = 9.8 - 0.005 V^2 m/s^2 for T0 = V0 = 0.0 and Tfinal = 10
Analytically, the integral of dv / ((1 - (0.005 / 9.8 ) V^2), which I don't want to have to look up (like surveying before the HP-35: eight place tables of Sin and Tan, and a Monroe hand-crank calculator). Great that there is a handy calculator function to solve those pesky DEs.

But the utility of the RKF solver is hampered save for Werner's guess that a << program object>> could be substituted for an 'algebraic object/expression', which is undocumented. I really do not think the firmware writers gave this much thought.

Thanks for the link to the pdf of HP-48 Insights.

[Giuseppe Donnini]

(06-15-2021 09:09 PM)rcthacker Wrote:  HP Solve doesn't solve DEs numerically.

I know, but that's beside the point. HP Solve is just one example of the general RPL feature we are talking about here, namely that for numeric evaluation any symbolic object may be replaced with a program object, provided the latter conforms to the input/output syntax of an algebraic expression. The same holds true for the plotter, for example, where any program equivalent to an expression (that is, which takes no arguments from the stack, and returns exactly one number) can be used instead of an expression.

(06-15-2021 09:09 PM)rcthacker Wrote:  Werner's guess that a << program object>> could be substituted for an 'algebraic object/expression', which is undocumented.

Presenting this general feature as a lucky “guess” on Werner's part verges on the absurd, when it has been clearly and officially documented since the times of the HP-28C. Indeed, in the HP-28C Reference Manual, which is available here, we may read the following:

Algebraic objects and programs have identical internal structures. Both types of procedures are sequences of objects that are processed sequentially when the procedures are evaluated. The algebraic 'X+Y' and the program « X Y + » are both stored as the same sequence (the RPN form). Algebraics are “marked” as algebraics so that they will be displayed as mathematical expressions and to indicate that they satisfy algebraic syntax rules. (p. 60)
[...]
If we take the above expression and rewrite it again by removing the parentheses, and placing the functions after their arguments, we obtain the RPN form of the expression. [...] This defines a program that has algebraic syntax, and is effectively equivalent to the corresponding algebraic object. (p. 62)

(06-15-2021 09:09 PM)rcthacker Wrote:  I really do not think the firmware writers gave this much thought.

And I really do not think that you do yourself a favor by calling people like Bill Wickes and Charlie Patton toughtless or sloppy.

In any event, here is some more background information on the subject by the firmware writers themselves (see “An Advanced Scientific Graphing Calculator” by Diana K. Byrne, Charles M. Patton, David Arnett, Ted W. Beers, and Paul J. McClellan, in: Hewlett-Packard Journal, Vol. 45 no. 4, August 1994, p. 21):

The HP 48G/GX contains two differential equation solvers and solution plotters. These solvers and solution plotters can be accessed via their input forms or invoked programmatically via commands. We provide a programmatic interface to the differential equation solvers and their subtasks so the user can use them with the calculator’s general solver feature to determine when a computed differential equation solution satisfies some condition, or to implement custom differential equation solvers from their subtasks.

In implementing the differential equation solution plots, one challenge was to identify and implement good solution methods. Another challenge was to merge this new plot type with the new 3D plot types described earlier and with the existing HP 48SX plot environment in a backward-compatible manner.

The HP 48G/GX specifically solves the initial value problem, consisting of finding the solution y(t) to the first-order equation y’(t) = f(t,y) with the initial condition y(t0) = y0. Here y'(t) denotes the first derivative of a scalar-valued or vector-valued solution y with respect to a scalar-valued parameter t. Higher-order differential equations can be expressed as a first-order system, so this problem is more general than it might at first appear.

Many solution methods have been developed over the years to solve the initial value problem. We decided to implement two methods, a Runge-Kutta-Fehlberg method for simplicity and speed of execution and a Rosenbrock method for reliability. The first method is easier to use, requiring less information from the user, but can fail on stiff problems. [*] The Rosenbrock method requires more information from the user, but can solve a wider selection of initial value problems. Both initial value problem solution methods require the user to provide the function f(t,y), the initial conditions, the final value of t, and an absolute error tolerance. The Rosenbrock method also requires the derivative of f(t,y) with respect to y (FYY) and the derivative of f(t,y) with respect to t (FYT).

All plot types use the contents of the variable EQ, typically to specify the function to be plotted. If the user selects the stiff (Rosenbrock) method the extra functions are passed to the solver by binding EQ to a list of functions f(t,y), FYY, and FYT. Otherwise, EQ is bound to the function f(t,y) needed by the Runge-Kutta-Fehlberg method.

Both methods solve the initial value problem by computing a series of solution steps from the initial conditions towards the final value, by default taking steps as large as possible subject to maintaining the specified error tolerance. The solution plotter plots the computed values and by default draws straight lines between the plotted points. However, although the computed steps may be accurate, the line segments drawn between the step endpoints may poorly represent the solution between those points. The plot parameter RES is used by many plot types to control the plot resolution. If RES is zero the initial value problem solution plotter imposes no additional limits on the step sizes. If RES is nonzero the plotter limits each step to have maximum size RES.

For the scalar-valued initial value problem it is typical to plot the computed solution y(t) on the vertical axis and the parameter t on the horizontal axis. However, in the vector-valued case the choice of what is to be plotted is not as clear. The user may wish a particular component of the computed solution plotted versus t or may wish two components plotted versus each other. The HP 48G/GX allows the user to specify the computed scalar solution, any component of the computed vector solution, or the parameter t to be plotted on either axis. This flexibility was introduced into the plot environment by expanding the AXES plot parameter. Previously, this parameter specified the coordinates of the axes origin. This parameter was expanded so that an optional form is a list specifying the origin and the horizontal and vertical plot components.

By judiciously expanding the meaning of the various plot parameters we were able to accommodate the differential equation solution plot type while maintaining backward compatibility with previous plot types.

[*] Stiff problems typically have solution components with large differences in time scale. More information is needed by a solver to compute a solution efficiently.

[jeejohn]

I have to agree with rcthacker about the "utility of the RKF solver".

I ran into this same issue some 20 odd years ago. The example in the HP 48G Series User's Guide showed the solution for a 2nd order differential equation with a Vector-Matrix equation. This works OK for Linear 2nd order differential equations or for a Linear System of 1st order differential equations. And even then within the SOLVE application - Differential Equations section.

Nowhere in any of the various books published for the 28 or 48s was an example given of solving a nonlinear 2nd order differential equation by way of using a program object.

On the other hand, I have not check to see how using Vector-Matrix equation on the 50G which allows for symbolic matrices might work.

John