1-AUG-2019 SolveSys - Equation Library and Nonlinear Systems Solver

[Han]

1-AUG-2019 (Ver. 1.10Beta)
CHANGELOG

• Fixed issues with fraction mark and separator characters. (Setting digits to 1234,56 using "," as the fraction mark would lead to unexpected behavior; should be fixed now).
• Preparing app for more inverse options (e.g. QR).
• Pressing [Shift][Plot] to access the library settings now has a new drop-down option to select which library to use. (Create new libraries by simply typing AFiles("My Library.lib"):={} and "My Library" will appear as an option.)
  

To update, first save your library files if you have modified any systems of equations. See the eqlib.pdf manual file for a program that exports library files into notes that can be edited on a computer. (The connectivity kit can also be used to view/edit notes.) It is always good practice to back up your libraries regularly.

Download links at bottom of post.

ABOUT THE APP



This app is intended to replicate and combine programs from the HP48 series of calculators for organizing equations and launching an equation solver. While there are a few user-created programs/apps that help manage equations, most of them have the formulas built-in (hard-coded into the source code of the program/app). This app attempts to allow users to build their own equations library. The user can then also launch a numerical solver. Using the built-in Solve app may be adequate for some users, but it lacks the ability to handle more general systems.

The current Solve app is quite deficient in that in requires the number of equations in a system be the same as the number of variables. Therefore, over-determined or under-determined systems cannot be solved numerically. Usually under-determine systems have an infinite number of results, and a numerical solver should ideally try to return a solution that is as close to the initial guess as possible. For over-determined systems, a solution usually does not exist. Most numerical solvers would return a "least-squares" solution.

This app implements the Newton-Raphson algorithm for solving a system of equations. In order to handle ill-conditioned, over-determined, under-determined, rank-deficient (etc.) Jacobian matrices, the SVD is used to determine the pseudo-inverse of the Jacobian (the pseudo-inverse coincides with the regular inverse of a matrix for non-singular matrices). Users have reported that the app can handle complex values just fine. Minimal testing seems to confirm this.

INSTALLATION AND USAGE

(Please refer to the reference manual included in the zip file)

REFERENCE MANUAL

The reference manual will likely see another update in which all the images are updated to nicer, crisp images created by \( \LaTeX\ \) and ps-tricks. For now, though, I will leave them there so that there is at least some frame of reference, and for old time's sake. The equations will also likely get a face lift as the variables are updated to follow typical mathematical conventions (e.g. subscripts). While the HP Prime currently does not support subscripted variables, the correspondence between the calculator's variables and those in the reference manual should be obvious (e.g. \( v_0 \) vs. \( \tt{v0} \) ).

If you notice any typos, please let me know.

OTHER NOTES:

There are two versions: a merged version (SVD merged into the source code for Equation Library) and the separated version. Does it make a difference which version I use? Yes, if you plan to work with under/over determined systems or systems whose Jacobian is either singular or near singular, it can make a major difference. The separated version (i.e. not merged) is slightly more accurate. As of the 2018 10 16 firmware release (press [Help] (Tree) (About) to see your firmware version), it appears that an issue with CAS corruption has been fixed. What this will likely mean is that the app will move toward a single release.

PRIOR CHANGES:
25-MARCH-2017 (Ver. 1.000)

• Major updates to reference manual
• Added support for animated diagrams (see Oscillations > Mass-Spring System)
• Library file format modified slightly (see ref. manual)
• Minor fix in library file for Conduction+Convection system (L1-L3 changed to l1-l3)
• Added option to revert to Solve app engine if anyone needs it
• Fixed deletion bug when library has more than 70 systems
• Systems may now be saved with diagram information
• Preliminary support for switching libraries (i.e. future update will allow users to have multiple libraries rather than saving everything into a single library)
• Information during solve process much more informative
• New icon (top half of old HP Solve application card)
  

08-MAR-2017 (Ver. 0.900)

• Added pictures; press [View] > View Picture
• Currently supports picture viewing; cannot save systems with a picture yet
• Fixed some typos in the reference manual
• This update modifies the Equation Library.lib file to add a 9th (picture data)
  

TODO

• Add a command to allow solving without using the app interface: ssSolve({ equations }, { vars }, { initial guesses })
• Add options to use other inverse techniques (e.g. QR, LU, or regular \( A^{-1} \)) to improve upon the solving speed
• Curve fitting like in SolveSys from the HP48/HP49/HP50 days
• Move all strings into a single list for localization.
• Implement numerical approximation of the Jacobian (right now, it is computed using the CAS which could potentially cause issues.
• Units
  

DOWNLOADS

Version 1.10 Beta: Attached in this recent post

Version 1.000 (same file name, just re-download to get the update)

  Equation Library Prime.zip (Size: 540.99 KB / Downloads: 2042)

[akmon]

Hello, downloaded beta version. Thank you for your efforts on try to have a solvesys like HP50, one of the most important programs, I think.

First, a strange issue has happend sometimes: after inserting the equations, I push [Num] then OK, and it does not start to solve but it returns to the equation menu. ¿? I did a warmstart and it didn´t happen again, at the moment.

The system of this thread gives different results from the Wolfram, I get A=1,27E-10 and B=1,088E21. Wolfram: A: -3,675E-10, B: 9,438€-11 Perhaps too low values for iteration accuracy?
http://www.hpmuseum.org/forum/thread-767...t=solvesys

Ok on this http://www.hpmuseum.org/forum/thread-614...t=solvesys

Ok on this http://www.hpmuseum.org/forum/thread-604...t=solvesys

Ok on this http://www.hpmuseum.org/forum/thread-238...t=solvesys

I´ll make more tests.

[Han]

(02-09-2017 08:42 PM)akmon Wrote:  Hello, downloaded beta version. Thank you for your efforts on try to have a solvesys like HP50, one of the most important programs, I think.

First, a strange issue has happend sometimes: after inserting the equations, I push [Num] then OK, and it does not start to solve but it returns to the equation menu. ¿? I did a warmstart and it didn´t happen again, at the moment.

The system of this thread gives different results from the Wolfram, I get A=1,27E-10 and B=1,088E21. Wolfram: A: -3,675E-10, B: 9,438€-11 Perhaps too low values for iteration accuracy?
http://www.hpmuseum.org/forum/thread-767...t=solvesys

If I recall correctly, that pair of values results in a nearly-singular Jacobian (likely singular due to floating point issues). But getting there was likely due to a few other issues. The SVD program had a few bugs in it (update again, even if you already did so recently). The other issue is due to CAS corruption, for which the current workaround is to use the restart command in the CAS view and rebuild the SVD program by opening and closing its source file.

I can still get both solutions to this system (in the 0.500 version) by setting |x| to 1e-20 and |F| to 1e-10. The |F| value I think should be self-explanatory but the |x| labeling for the other tolerance may be a bit misleading. It is really a measure of the change in the x_i's relative to the norm of the previous solution. So if u want 4 digits of accuracy and your solution is in the order of 1e-10, then you should try setting |x| to 1e-16. May I can change this later on to "number of digits" rather than a "tolerance" value to reflect the true behavior of that value. Also, |F| should be set to be a 1e-10 (or smaller) since 1.27e-10 would be pretty close to the error itself (similar reasoning). When using (1,1) and (1,2) for initial values, I get the two separate solutions.

Lastly, keep in mind that the (now prior) version was just Newton's method and is highly likely to wander off into "la-la" land. In your case, you probably landed in a spot where the Jacobian is nearly singular (or in fact singular but not shown to be so due to typical floating point issues). You can check this by using svd2(ssJ) -- make sure to do this in the CAS view. The "middle" vector lists the singular values. If any of them are close to zero, then the Jacobian is (likely) rank-deficient.

The newer version has backtracking line search.

I still need to tweak the exit conclusions (i.e. how to classify supposed solutions).

I just updated the program; I appreciate your continued efforts to test the solver.

EDIT: I should mention that I am able to get it to run on the real calculator after having tweaked a few spots where iferr blocks were used. It would seem that iferr and return calls are currently misbehaving so that repeated instances would lead to a warmstart. I worked around that issue by modifying the code to not rely on iferr.

[akmon]

Thank you for the new release, I have also updated SVD. I have edited the solver settings with your advice, edited the solver settings, previous values [1,1], and after a few seconds I get A=3,36E-5 ; B= 1,776€-5, (100 iterations). After the results I started to solve again with these new initial values, and after a few seconds the result was correct. BUT, I don´t know why I cannot exit from the matrix environment. I can go to the solver menu, but no more, no Home, no CAS, strange, only happened once, tried to do the same again and no problem this time.

Edited solver parameters to 200 iterations and the result was correct.

and if I assign [1,2] it shows "linesearch failed f> fold): yes lam 0,371136323793
IncrX=7,45499482206

If I use [-1,1] or [0,0] or [1,0] and press OK, it returms to equation menu. ¿?

By the way, the clocksand didn´t appear while solving with [1,1], I thought it hanged, but finished the calculus.
And it seems that sometimes SVD does not exist. I edit the program or debug, and it works again.

[Han]

(02-10-2017 12:18 AM)akmon Wrote:  Thank you for the new release, I have also updated SVD. I have edited the solver settings with your advice, edited the solver settings, previous values [1,1], and after a few seconds I get A=3,36E-5 ; B= 1,776€-5, (100 iterations). After the results I started to solve again with these new initial values, and after a few seconds the result was correct. BUT, I don´t know why I cannot exit from the matrix environment. I can go to the solver menu, but no more, no Home, no CAS, strange, only happened once, tried to do the same again and no problem this time.

Edited solver parameters to 200 iterations and the result was correct.

and if I assign [1,2] it shows "linesearch failed f> fold): yes lam 0,371136323793
IncrX=7,45499482206

If I use [-1,1] or [0,0] or [1,0] and press OK, it returms to equation menu. ¿?

By the way, the clocksand didn´t appear while solving with [1,1], I thought it hanged, but finished the calculus.
And it seems that sometimes SVD does not exist. I edit the program or debug, and it works again.

It looks like you have memory corruption (referring to the bold text).

The solving process is a bit slower in this version due to the backtracking. Before, every iteration step was a full Newton step in the opposite direction of the gradient. This could mean fewer iterations, and in other cases it could mean the Newton steps wander around in an unwanted direction. The backtracking line search technique takes smaller steps, but ensures that each step is a descent step toward either a 0 or a minimum. This is the tradeoff; perhaps I can make it optional to use line searches. As for the lockup issues:

Try the following:

1) type restart into the CAS view
2) go into the memory browser, and look in the CAS Vars submenu; delete anything that is still left over in this submenu -- make sure this is completely clean
3) turn the calculator off, and then back on (this presumably forces the memory to save its state into flash)
4) warmstart using [On][Symb]

When the calculator warmstarts, it should rebuild all stored programs. If you see a message about installing SVD then open and immediately close the source file for SVD.

Then try your solve.

[akmon]

Did all the steps you suggested, but the results are exactly the same as the post I worte yesterday, the same behaviour (tests on my physical calculator). If you need any parameter about my settings, please let me know.

[Han]

It looks like our equations are slightly different. If your equation is \( \left( \frac{A+B}{A} \right) ^2 = 8.377 \) then a guess of A=1 and B=2 would lead to a solution whereas A=2 and B=1 would lead to a failed line search.

If your equation is \( \left( \frac{A+B}{B} \right) ^2 = 8.377 \) -- notice the denominator -- then a guess of A=2 and B=1 would lead to a solution; A=1 and B=2 would lead to a failed line search (symmetry in the variables compared to the first case)

For the rest of this post, let us assume you used the first form of (with A in the denominator). Assuming that our application is bug-free (doubtful at the moment), then A=2 and B=1 is a point for which a line search is going to fail. To see how this might come to be, let me explain the algorithm a bit.

In the recent versions, we are additionally solving for the minimums of another function f = 1/2 (F.F). Notice that the zeros of F are guaranteed minimums of f (since it is a non-negative function). However, not all minimums of f are zeros of F -- no worries, we simply check whether \( \| F \| \) is sufficiently small enough). We are still using Newton's method, wherein we compute \( \Delta x = -J^{-1} \cdot F \) and obtain \( x_\text{new} = x_\text{old} + \Delta x \), however we scale the step \( \Delta x \) by \( \lambda \) so that we guarantee a decrease in \( \| F \| \). (J here represents the Jacobian of F). This is the line search portion of the algorithm.

You can read about line searches here:

http://reference.wolfram.com/language/tu...thods.html

Basically, the line search tries to find a suitable \( \lambda \) that satisfies the strong Wolfe conditions (so that we have convergence). (In the link above, they call the scalar factor \( \alpha^* \) ). From the link:

Quote:Typically, it only takes a few iterations to find a point satisfying the conditions. However, the line search is not always able to find a point that satisfies the conditions. Usually this is because there is insufficient precision to compute the points closely enough to satisfy the conditions, but it can also be caused by functions that are not completely smooth or vary extremely slowly in the neighborhood of a minimum.

So in the end, A=2 and B=1 is just not a good guess (again, assuming A is the variable in the denominator of your equation involving 8.377).

[akmon]

Understood. I have always made the examples using B as denominator, so if I insert A=1 and B=2 then an error message appears as yours. But with A=2 and B=1 I find the solution after 180 iterations
Please use the equations with B denominator and Guess with [0,0] or [1,0], it returns to equation menu again. As I have now in Verbose mode, I can read (slightly) that before going to equation menu, the program ends in creating the F vector.

On virtual prime is even worse. It crashes with those guess values.

[Han]

(02-10-2017 07:13 PM)akmon Wrote:  On virtual prime is even worse. It crashes with those guess values.

Wow, that is a huge crash!

[Han]

I think to get around all these firmware issues, I will be

1) converting the svd code from CAS form into non-CAS form
2) then embedding it within the equation library / solver app

This will mean no more corruption of the svd2 program (although a few variables must unfortunately remain on the cas side) and not having to update two separate programs. There is a tiny downside (possibly) since the home programs use BCD representation of floating points and the cas programs use binary representation. I am not sure how many bits are used on either side, but it could mean a tiny bit of accuracy loss switching out of cas form.

[Han]

See if this version works any better. It will still misbehave in the event of something like you described above (and crash the VC), but this version no longer makes use of CAS programs to do the svd calculation. All of that is now built-in and part of the app as a whole. So you will need to delete your existing copy of svd2 as they are now part of the app .

For good measure, clean out your CAS variables to ensure none of the old stuff is still lying around, possibly causing collisions with their duplicate versions inside the app's code.


  eqlibsvd_10feb2017.txt (Size: 62.59 KB / Downloads: 62)

Anyone interested in creating a nice icon for this app? Or perhaps creating the data file that would correspond to the equation library from the HP48 series?

EDIT: The non-CAS version of svd2 appears to be slightly less accurate than the CAS version (off at about the 12th digit) (I'm guessing this is due to the floating point format)

[akmon]

Now that SVD is included in the core of the app, it´s not only more conformtable for updating, but it avoids errors from "SVD not found". Well done. Buy the way, remember to update your first post.

I haved noticed that slightly less accuracy in our famous system (A= -3,4985E-10, B=8,98E-11), I have also discovered that if you are very accurate on your guess you´ll reach the exact solution, but it is too much to suppose and this system of equations are very unlikely to appear in general terms. What still concern me is what causes that fatal crash on VC just only choosing from one guess to other.

By the way, I´ve tried to solve a complex system. Always get bad guesses. How do I solve this? (2+i)*X+2*Y=1+7*i ; (1-i)*X+i*Y=0

About icons, three examples. The first, represents the graphical purpose of this app. The second, a wink to the firs Solvesys, and the last, another wink to the first EQlib from HP48.

[Han]

(02-12-2017 11:19 PM)akmon Wrote:  Now that SVD is included in the core of the app, it´s not only more conformtable for updating, but it avoids errors from "SVD not found". Well done. Buy the way, remember to update your first post.

I haved noticed that slightly less accuracy in our famous system (A= -3,4985E-10, B=8,98E-11), I have also discorvered that if you are very accurate on your guess you´ll reach the exact solution, but it is too much to suppose and this system of equations are very unlikely to appear in general terms. What still concern me is what causes that fatal crash on VC just only choosing from one guess to other.

This is exactly what I was referring to regarding the loss in accuracy. Using the CAS version, you get just a slight advantage in accuracy. (see http://www.hpmuseum.org/forum/thread-2953.html)

But you are right, such a system tends to be atypical due to the singularity at B=0 (assuming (((A+B)/B)^2\) ), or A=-B

Right now, you should only need to update the SVD program if there's a bug fix. However, you will still have to deal with the CAS corruption (esp. on the emulator). On the hardware, it should not occur unless you either 1) recompiled a program or 2) experience a warmstart. My posting the code as source code probably doesn't help with (1).

The crash is due to the firmware's inability to handle undefined values in the Home view. The app is programmed as a "non-CAS" program. However, it still makes a few calls into the CAS. I hope that they will fix this in the next firmware.

Quote:About icons, three examples. The first, represents the graphical purpose of this app. The second, a wink to the firs Solvesys, and the last, another wink to the first EQlib from HP48.

Ahh, I do like the nod to either SolveSys or EQLIB. And compsystems' idea about using the old solve card as an icon is also a nice idea.

[Han]

Equation Library data (partially completed). With the data below, you can get an idea of what the final version might look like. Please note that it is merely an example of a .lib file and is currently incomplete, or even wrong in some parts as I just recently learned.

If you would like a data file to test, enter the app and copy the list of lists below (easiest if done on emulator), and type:

AFiles("Equation Library.lib"):=<insert copied list of lists below here>

Code:


{
{
  "Elastic Buckling",
  {
    "Pcr=π^2*E*A/((K*L)/r)^2",
    "Pcr=π^2*E*I/(K*L)^2",
    "σcr=Pcr/A",
    "r=√(I/A)"
  },
  { 1,1,1,1 },
  {
    "L", "K", "r", "E", "I", "A", "Pcr", "σcr"
  },
  {
    "Length of column",
    "Effective length factor",
    "Radius of gyration",
    "Modulus of elasticity",
    "Moment of inertia",
    "Cross-sectional area",
    "Critical load",
    "Critical stress"
  },
  { 0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Eccentric Columns",
  {
    "σmax=P/A*(1+ε*c/r^2*(1/COS(K*L/(2*r)*√(P/(E*A)))))",
    "r=√(I/A)"
  },
  { 1,1 },
  {
    "L", "K", "r", "E", "I", "A",
    "ε", "c", "σmax", "P"
  },
  {
    "Length of column",
    "Effective length factor",
    "Radius of gyration",
    "Modulus of elasticity",
    "Moment of Inertia",
    "Cross-sectional area",
    "Eccentricity of load",
    "Distance to edge fiber",
    "Maximum stress",
    "Load"
  },
  { 0,0,0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Simple Deflection",
  { 
    "y=P*(L-a)*x/(6*L*E*I)*(x^2+(L-a)^2-L^2)-M*x/(E*I)*(c-x^2/(6*L)-L/3-c^2/(2*L))-w*x/(24*E*I)*(L^3+x^2*(x-2*L))"
  },
  { 1 },
  {
    "y", "P", "L", "a", "x", "E", "I", "M", "c", "w"
  },
  {
    "Deflection at x",
    "Point load",
    "Length of beam",
    "Distance to point load",
    "Distance along beam",
    "Modulus of elasticity",
    "Moment of inertia",
    "Applied moment",
    "Dist. to applied moment",
    "Distributed load"
  },
  { 0,0,0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Simple Slope",
  {
    "Θ=P*(L-a)/(6*L*E*I)*(3*x^2+(L-a)^2-L^2)-M/(E*I)*(c-x^2/(2*L)-L/3-c^2/(2*L))-w/(24*E*I)*(L^3+x^2*(4*x-6*L))"
  },
  { 1 },
  {
    "Θ", "P", "L", "a", "x", "E", "I", "M", "c", "w"
  },
  {
    "Slope x",
    "Point load",
    "Length of beam",
    "Distance to point load",
    "Distance along beam",
    "Modulus of elasticity",
    "Moment of inertia",
    "Applied moment",
    "Dist. to applied moment",
    "Distributed load"
  },
  { 0,0,0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Simple Moment",
  {
    "Mx=P*(L-a)*x/L+M*x/L+w*x/2*(L-x)"
  },
  { 1 },
  {
    "Mx", "P", "L", "a", "x", "M", "w"
  },
  {
    "Internal moment at x",
    "Point load",
    "Length of beam",
    "Cross-sectional area",
    "Distance along beam",
    "Applied moment",
    "Distributed load"
  },
  { 0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Simple Shear",
  {
    "V=P*(L-a)/L+M/L+w/2*(L-2*x)"
  },
  { 1 },
  {
    "V", "P", "L", "a", "M", "w", "x"
  },
  {
    "Shear force at x",
    "Point load",
    "Length of beam",
    "Cross-sectional area",
    "Applied moment",
    "Distributed load",
    "Distance along beam"
  },
  { 0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Cantilever Deflection",
  {
    "y=P*x^2/(6*E*I)*(x-3*a)+M*x^2/(2*E*I)-w*x^2/(24*E*I)*(6*L^2-4*L*x+x^2)"
  },
  { 1 },
  {
    "y", "P", "L", "a", "x", "E", "I", "M", "w"
  },
  {
    "Deflection at x",
    "Point load",
    "Length of beam",
    "Distance to point load",
    "Distance along beam",
    "Modulus of elasticity",
    "Moment of inertia",
    "Applied moment",
    "Distributed load"
  },
  { 0,0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Cantilever Slope",
  {
    "Θ=P*x/(2*E*I)*(x-2*a)+M*x/(E*I)-w*x/(6*E*I)*(3*L^2-3*L*x+x^2)"
  },
  { 1 },
  {
    "Θ", "P", "L", "a", "x", "E", "I", "M", "w"
  },
  {
    "Slope x",
    "Point load",
    "Length of beam",
    "Distance to point load",
    "Distance along beam",
    "Modulus of elasticity",
    "Moment of inertia",
    "Applied moment",
    "Distributed load"
  },
  { 0,0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Cantilever Moment",
  {
    "Mx=P*(x-a)+M-w/2*(L^2-2*L*x+x^2)"
  },
  { 1 },
  {
    "Mx", "P", "L", "a", "x", "M", "w"
  },
  {
    "Internal moment at x",
    "Point load",
    "Length of beam",
    "Cross-sectional area",
    "Distance along beam",
    "Applied moment",
    "Distributed load"
  },
  { 0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Cantilever Shear",
  {
    "V=P+w*(L-x)"
  },
  { 1 },
  {
    "V", "P", "L", "x", "w"
  },
  {
    "Shear force at x",
    "Point load",
    "Length of beam",
    "Distance along beam",
    "Distributed load"
  },
  { 0,0,0,0,0 },
  { 0,0,0,0,0 },
  "Columns and Beams"
},

{
  "Coulomb's Law",
  {
    "F=1/(4*π*ε0*εr)*(q1*q2/r^2)"
  },
  { 1 },
  {
    "q1", "q2", "r", "εr", "F", "ε0"
  },
  {
    "Point 1 charge",
    "Point 2 charge",
    "Charge distance",
    "Relative permittivity",
    "Force between charges",
    "(Permittivity Constant)"
  },
  { 0,0,0,0,0,8.85418781762e-12 },
  { 0,0,0,0,0,1 },
  "Electricity"
},

{
  "Ohm's Law and Power",
  {
    "V=I*R",
    "P=V*I",
    "P=I^2*R",
    "P=V^2/R"
  },
  { 1,1,1,1 },
  {
    "V", "I", "R", "P"
  },
  {
    "Voltage",
    "Current",
    "Resistance",
    "Power"
  },
  { 0,0,0,0 },
  { 0,0,0,0 },
  "Electricity"
},

{
  "Voltage Divider",
  {
    "V1=V*(R1/(R1+R2))"
  },
  { 1 },
  {
    "V1", "V", "R1", "R2"
  },
  {
    "Voltage across R1",
    "Total voltage",
    "Resistance 1",
    "Resistance 2"
  },
  { 0,0,0,0 },
  { 0,0,0,0 },
  "Electricity"
},

{
  "Current Divider",
  {
    "I1=I*(R2/(R1+R2))"
  },
  { 1 },
  {
    "I1", "I", "R1", "R2"
  },
  {
    "Current in R1",
    "Total current",
    "Resistance 1",
    "Resistance 2"
  },
  { 0,0,0,0 },
  { 0,0,0,0 },
  "Electricity"
},

{
  "Wire Resistance",
  {
    "R=ρ*L/A"
  },
  { 1 },
  {
    "R", "ρ", "L", "A"
  },
  {
    "Resistance",
    "Resistivity",
    "Length of wire",
    "Cross-sectional area"
  },
  { 0,0,0,0 },
  { 0,0,0,0 },
  "Electricity"
},

{
  "Series/Parallel Resistance",
  {
    "Rs=R1+R2",
    "1/Rp=1/R1+1/R2"
  },
  { 1,1 },
  {
    "Rs", "Rp", "R1", "R2"
  },
  {
    "Series resistance",
    "Parallel resistance",
    "Resistance 1",
    "Resistance 2"
  },
  { 0,0,0,0 },
  { 0,0,0,0 },
  "Electricity"
},

{
  "Series/Parallel Capacitance",
  {
    "Cp=C1+C2",
    "1/Cs=1/C1+1/C2"
  },
  { 1,1 },
  {
    "Cs", "Cp", "C1", "C2"
  },
  {
    "Series capacitance",
    "Parallel capacitance",
    "Capacitance 1",
    "Capacitance 2"
  },
  { 0,0,0,0 },
  { 0,0,0,0 },
  "Electricity"
},

{
  "Series/Parallel Inductance",
  {
    "Ls=l1+l2",
    "1/Lp=1/l1+1/l2"
  },
  { 1,1 },
  {
    "Ls", "Lp", "l1", "l2"
  },
  {
    "Series inductance",
    "Parallel inductance",
    "Inductance 1",
    "Inductance 2"
  },
  { 0,0,0,0 },
  { 0,0,0,0 },
  "Electricity"
},

{
  "Capacitive Energy",
  {
    "E=C*V^2/2"
  },
  { 1 },
  {
    "E", "C", "V"
  },
  {
    "Energy",
    "Capacitance",
    "Voltage"
  },
  { 0,0,0 },
  { 0,0,0 },
  "Electricity"
},

{
  "Inductive Energy",
  {
    "E=L*I^2/2"
  },
  { 1 },
  {
    "E", "L", "I"
  },
  {
    "Energy",
    "Inductance",
    "Current"
  },
  { 0,0,0 },
  { 0,0,0 },
  "Electricity"
},

{
  "RLC Current Delay",
  {
    "TAN(φs)=(XL-XC)/R",
    "TAN(φp)=(1/XC-1/XL)/(1/R)",
    "XC=1/(ω*C)",
    "XL=ω*L",
    "ω=2*π*f"
  },
  { 1,1,1,1,1 },
  {
    "ω", "f", "C", "L", "XC", "XL", "φs", "φp"
  },
  {
    "Angular frequency",
    "Frequency",
    "Capacitance",
    "Inductance",
    "Reactance of capacitor",
    "Reactance of inductor",
    "Series phase angle",
    "Parallel phase angle"
  },
  { 0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0 },
  "Electricity"
},

{
  "DC Capacitor Current",
  {
    "I=C*(ΔV/Δt)",
    "ΔV=Vf-Vi",
    "Δt=tf-ti"
  },
  { 1,1,1 },
  {
    "I", "C", "ΔV", "Vf", "Vi", "Δt", "tf", "ti"
  },
  {
    "Current",
    "Capacitance",
    "Voltage change",
    "Final voltage",
    "Initial voltage",
    "Time change",
    "Final time",
    "Initial time"
  },
  { 0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0 },
  "Electricity"
},

{
  "Capacitor Charge",
  { "q=C*V" },
  { 1 },
  { "q", "C", "V" },
  {
    "Charge",
    "Capacitance",
    "Voltage"
  },
  { 0,0,0 },
  { 0,0,0 },
  "Electricity"
},

{
  "DC Inductor Voltage",
  {
    "V=-L*(ΔI/Δt)",
    "ΔI=I1-I0",
    "Δt=t1-t0"
  },
  { 1,1,1 },
  {
    "V", "L", "ΔI", "I1", "I0", "Δt", "t1", "t0"
  },
  {
    "Voltage",
    "Inductance",
    "Current change",
    "Final current",
    "Initial current",
    "Time change",
    "Final time",
    "Initial time"
  },
  { 0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0 },
  "Electricity"
},

{
  "RC Transient",
  { "V=Vf-(Vf-Vi)*e^((-t)/(R*C))" },
  { 1 },
  {
    "V", "Vf", "Vi", "t", "R", "C"
  },
  {
    "Capacitor voltage",
    "Final voltage",
    "Initial voltage",
    "Time",
    "Resistance",
    "Capacitance"
  },
  { 0,0,0,0,0,0 },
  { 0,0,0,0,0,0 },
  "Electricity"
},

{
  "RL Transient",
  { "I=1/R*(Vf-(Vf-Vi)*e^(-t*R/L))" },
  { 1 },
  {
    "I", "R", "Vf", "Vi", "t", "L"
  },
  {
    "Inductor Current",
    "Resistance", 
    "Final voltage",
    "Initial voltage",
    "Time",
    "Inductance"
  },
  { 0,0,0,0,0,0 },
  { 0,0,0,0,0,0 },
  "Electricity"
},

{
  "Resonant Frequency",
  {
    "ω0=1/√(L*C)",
    "Qs=1/R*√(L/C)",
    "Qp=R*√(C/L)",
    "ω0=2*π*f0"
  },
  { 1,1,1,1 },
  {
    "ω0", "L", "C", "Qs", "Qp", "R", "f0"
  },
  {
    "Resonant angular frequency",
    "Inductance",
    "Capacitance",
    "Series quality factor",
    "Parallel quality factor",
    "Resistance",
    "Resonant frequency"
  },
  { 0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0 },
  "Electricity"
},

{
  "Plate Capacitor",
  { "C=(ε0*εr*A)/d" },
  { 1 },
  {
    "C", "εr", "A", "d", "ε0"
  },
  {
    "Capacitance",
    "Relative permittivity",
    "Plate area",
    "Plate separation",
    "(Permittivity Constant)"
  },
  { 0,0,0,0,8.85418781762e-12 },
  { 0,0,0,0,1 },
  "Electricity"
},

{
  "Cylindrical Capacitance",
  { "C=2*π*ε0*εr*L/LN(ro/ri)" },
  { 1 },
  {
    "C", "εr", "L", "ro", "ri", "ε0"
  },
  {
    "Capacitance",
    "Relative permittivity",
    "Length of cylinder",
    "Outside radius",
    "Inside radius",
    "(Permittivity Constant)"
  },
  { 0,0,0,0,0,8.85418781762e-12 },
  { 0,0,0,0,0,1 },
  "Electricity"
},

{
  "Solenoidal Inductance",
  { "L=μ0*μr*n^2*A*h" },
  { 1 },
  {
    "L", "μr", "n", "A", "h", "μ0"
  },
  {
    "Inductance",
    "Relative permeability",
    "Turns per unit length",
    "Cross-sectional area",
    "Length of solenoid",
    "(Permeability Constant)"
  },
  { 0,0,0,0,0,1.25663706144e−6},
  { 0,0,0,0,0,1 },
  "Electricity"
},  

{
  "Toroid Inductance",
  { "L=(μ0*μr*N^2*h)/(2*π)*LN(ro/ri)" },
  { 1 },
  {
    "L", "μr", "N", "h", "ro", "ri", "μ0"
  },
  {
    "Inductance",
    "Relative permeability",
    "Total number of turns",
    "Toroid thickness",
    "Outer radius",
    "Inner radius",
    "(Permeability Constant)"
  },
  { 0,0,0,0,0,0,1.25663706144e−6 },
  { 0,0,0,0,0,0,1 },
  "Electricity"
},

{
  "Sinusoidal Voltage",
  {
    "V=Vmax*SIN(ω*t+φ)",
    "ω=2*π*f"
  },
  { 1,1 },
  {
    "V", "Vmax", "ω", "φ", "t", "f"
  },
  {
    "Voltage",
    "Maximum voltage",
    "Angular frequency",
    "Phase angle",
    "Time",
    "Frequency"
  },
  { 0,0,0,0,0,0 },
  { 0,0,0,0,0,0 },
  "Electricity"
},

{
  "Sinusoidal Current",
  {
    "I=Imax*SIN(ω*t+φ)",
    "ω=2*π*f"
  },
  { 1,1 },
  {
    "I", "Imax", "ω", "φ", "t", "f"
  },
  {
    "Current",
    "Maximum current",
    "Angular frequency",
    "Phase angle",
    "Time",
    "Frequency"
  },
  { 0,0,0,0,0,0 },
  { 0,0,0,0,0,0 },
  "Electricity"
},

{
  "Pressure at Depth",
  { "P=P0+ρ*g*h" },
  { 1 },
  {
    "P", "P0", "ρ", "h", "g"
  },
  {
    "Pressure",
    "Reference pressure",
    "Density",
    "Depth",
    "(Acceleration Constant)"
  },
  { 0,0,0,0,9.80665 },
  { 0,0,0,0,1 },
  "Fluids"
},

{
  "Bernoulli Equation",
  {
    "ΔP/ρ+(v2^2-v1^2)/2+g*Δy=0",
    "ΔP/ρ+v2^2*(1-(A2/A1)^2)/2+g*Δy=0",
    "ΔP/ρ+v1^2*((A1/A2)^2-1)/2+g*Δy=0",
    "ΔP=P2-P1",
    "Δy=y2-y1",
    "M=ρ*Q",
    "Q=A2*v2",
    "Q=A1*v1",
    "A1=(π*D1^2)/4",
    "A2=(π*D2^2)/4"
  },
  { 1,1,1,1,1,1,1,1,1,1 },
  {
    "ΔP", "ρ", "v1", "v2", "Δy", "y1", "y2", "A1", "A2", "Q",
    "P1", "P2", "D1", "D2", "M", "g"
  },
  {
    "Pressure change",
    "Density",
    "Initial velocity",
    "Final velocity",
    "Height change",
    "Initial height",
    "Final height",
    "Initial area",
    "Final area",
    "Volume flow rate",
    "Initial pressure",
    "Final pressure",
    "Initial diameter",
    "Final diameter",
    "Mass flow rate",
    "(Acceleration Constant)"
  },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9.80665 },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1 },
  "Fluids"
},

{
  "Flow with Losses",
  {
    "M*(ΔP/ρ+(v2^2-v1^2)/2+g*Δy+hL)=W",
    "M*(ΔP/ρ+v2^2*(1-(A2/A1)^2)/2+g*Δy+hL)=W",
    "M*(ΔP/ρ+v1^2*((A1/A2)^2-1)/2+g*Δy+hL)=W",
    "ΔP=P2-P1",
    "Δy=y2-y1",
    "M=ρ*Q",
    "Q=A2*v2",
    "Q=A1*v1",
    "A1=(π*D1^2)/4",
    "A2=(π*D2^2)/4"
  },
  { 1,1,1,1,1,1,1,1,1,1 },
  {
    "ΔP", "ρ", "v1", "v2", "Δy", "y1", "y2", "A1", "A2", "Q",
    "P1", "P2", "D1", "D2", "M", "W", "hL", "g"
  },
  {
    "Pressure change",
    "Density",
    "Initial velocity",
    "Final velocity",
    "Height change",
    "Initial height",
    "Final height",
    "Initial area",
    "Final area",
    "Volume flow rate",
    "Initial pressure",
    "Final pressure",
    "Initial diameter",
    "Final diameter",
    "Mass flow rate",
    "Power input",
    "Head loss",
    "(Acceleration Constant)"
  },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9.80665 },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1 },
  "Fluids"
},

{
  "Flow in Full Pipes",
  {
    "ρ*(π*D^2/4)*vavg*(ΔP/ρ+g*Δy+vavg^2*(2*f*(L/D)+ΣK/2))=W",
    "ΔP=P2-P1",
    "Δy=y2-y1",
    "M=ρ*Q",
    "Q=A*vavg",
    "A=(π*D^2)/4",
    "RN=D*vavg*ρ/μ",
    "n=μ/ρ"
  },
  { 1,1,1,1,1,1,1,1 },
  {
    "ρ", "ΔP", "P1", "P2", "Δy", "y1", "y2", "f", "L", "D",
    "vavg", "W", "ΣK", "M", "RN", "Q", "A", "μ", "n", "g"
  },
  {
    "Density",
    "Pressure change",
    "Initial Pressure",
    "Final Pressure",
    "Height change",
    "Initial height",
    "Final height",
    "f ???",
    "Length",
    "Diameter",
    "Average velocity",
    "Power input",
    "Total fitting coefficients",
    "Mass flow rate",
    "Reynolds number",
    "Volume flow rate",
    "Cross-sectional area",
    "Dynamic viscosity",
    "Kinematic viscosity",
    "(Accelerations Constant)"
  },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9.80665 },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1 },
  "Fluids"
},

{
  "Linear Mechanics",
  {
    "F=m*a",
    "Ki=1/2*m*vi^2",
    "Kf=1/2*m*vf^2",
    "W=F*x",
    "W=Kf-Ki",
    "P=F*v",
    "Pavg=W/t",
    "vf=vi+a*t"
  },
  { 1,1,1,1,1,1,1,1 },
  {
    "F", "m", "a", "Ki", "Kf", "vi", "vf", "W", "x", "P", 
    "Pavg", "v", "t"
  },
  {
    "Force",
    "Mass",
    "Acceleration",
    "Initial kinetic energy",
    "Final kinetic energy",
    "Initial velocity",
    "Final velocity",
    "Work",
    "Distance",
    "Instantaneous power",
    "Average power",
    "Velocity",
    "Time"
  },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0 },
  "Forces and Energy"
},

{
  "Angular Mechanics",
  {
    "τ=I*α",
    "Ki=1/2*I*ωi^2",
    "Kf=1/2*I*ωf^2",
    "W=τ*θ",
    "W=Kf-Ki",
    "P=τ*ω",
    "Pavg=W/t",
    "ωf=ωi+α*t",
    "ta=α*r",
    "ω=2*π*N",
    "ωi=2*π*Ni",
    "ωf=2*π*Nf"
  },
  { 1,1,1,1,1,1,1,1,1,1,1,1 },
  {
    "τ", "I", "α", "Ki", "ωi", "Kf", "ωf", "W", "θ", "P",
    "Pavg", "t", "N", "Ni", "Nf", "r", "ω", "ta"
  },
  {
    "Torque",
    "Moment of inertia",
    "Angular acceleration",
    "Initial kinetic energy",
    "Initial angular velocity",
    "Final kinetic energy",
    "Final angular velocity",
    "Work",
    "Angular displacement",
    "Instantaneous power",
    "Average power",
    "Time",
    "Rotational speed",
    "Initial rotational speed",
    "Final rotational speed",
    "Radius",
    "Angular velocity",
    "Tangential acceleration"
  },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
  "Forces and Energy"
},

{
  "Centripetal Force",
  {
    "F=m*ω^2*r",
    "ω=v/r",
    "ar=v^2/r",
    "ω=2*π*N"
  },
  { 1,1,1,1 },
  {
    "F", "m", "ω", "r", "v", "ar", "N"
  },
  {
    "Force",
    "Mass",
    "Angular velocity",
    "Radius",
    "Velocity",
    "Centripetal acceleration",
    "Rotational speed"
  },
  { 0,0,0,0,0,0,0 },
  { 0,0,0,0,0,0,0 },
  "Forces and Energy"
},

{
  "Hooke's Law",
  {
    "F=-k*x",
    "W=-1/2*k*x^2"
  },
  { 1,1 },
  {
    "F", "k", "W", "x"
  },
  {
    "Force",
    "Spring constant",
    "Work",
    "Displacement"
  },
  { 0,0,0,0 },
  { 0,0,0,0 },
  "Forces and Energy"
},

{
  "1D Elastic Collisions",
  {
    "v1f=(m1-m2)/(m1+m2)*v1i",
    "v2f=(2*m1)/(m1+m2)*v1i"
  },
  { 1,1 },
  {
    "v1f", "v1i", "v2f", "v2i", "m1", "m2"
  },
  {
    "Object 1 final velocity",
    "Object 1 initial velocity",
    "Object 2 final velocity",
    "Object 2 initial velocity",
    "Object 1 mass",
    "Object 2 mass"
  },
  { 0,0,0,0,0,0 },
  { 0,0,0,0,0,0 },
  "Forces and Energy"
},

{
  "Drag Force",
  { "F=Cd*((ρ*v^2)/2)*A" }, 
  { 1 },
  {
    "F", "Cd", "ρ", "v", "A"
  },
  {
    "Drag force",
    "Drag coefficient",
    "Fluid density",
    "Velocity",
    "Projected area"
  },
  { 0,0,0,0,0 },
  { 0,0,0,0,0 },
  "Forces and Energy"
},

{
  "Law of Gravitation",
  { "F=G*((m1*m2)/r^2)" },
  { 1 },
  {
    "F", "m1", "m2", "r", "G"
  },
  {
    "Force",
    "Mass 1",
    "Mass 2",
    "Separation distance",
    "Attractive force",
    "(Gravitational Constant)"
  },
  { 0,0,0,0,6.67384e-11 },
  { 0,0,0,0,1 },
  "Forces and Energy"
},

{
  "Mass-Energy Relation",
  { "E=m*c^2" },
  { 1 },
  { "E", "m", "c" },
  {
    "Energy",
    "Mass",
    "(Speed of Light Constant)"
  },
  { 0,0,299792458 },
  { 0,0,1 },
  "Forces and Energy"
}
}

[Han]

Utility to convert old library data to newer format:

Code:

EXPORT UpdateLibDat()
BEGIN
  local lib:=AFiles("Equation Library.lib");
  local n:=size(lib);
  local j,sys;
  if n then
    for j from 1 to n do
      sys:=lib(j);
      if size(sys)==7 then
        sys(0):="New Category";
      end;
      lib(j):=sys;
    end;
    AFiles("Equation Library.lib"):=lib;
  end;
END;

[Han]

Here is another utility that will take your library data and convert it into a note, which you can then export into a text file. It assumes your library data has be converted up to the version used in 0.651

Code:

EXPORT LibToNote()
BEGIN
  local note:="EqLib.txt";
  local lib:=AFiles("Equation Library.lib");
  local n:=size(lib);
  local sys, m, j, k, r;
  local indent;

  Notes(note):="{\n";
  for j from 1 to n do
    Notes(note):=Notes(note) + "  {\n";
    sys:=lib(j);
    for k from 1 to 8 do
      Notes(note):=Notes(note) + "    ";
      if (k == 2) or (k == 5) then
        Notes(note):=Notes(note) + "{\n";
        for r from 1 to (size(sys(k))-1) do
          Notes(note):=Notes(note) + "      ";
          Notes(note):=Notes(note) + string(sys(k,r)) + ",\n";
        end;
        Notes(note):=Notes(note) + "      ";
        Notes(note):=Notes(note) + string(sys(k,r)) + "\n    }";
      else
        Notes(note):=Notes(note) + string(sys(k));
      end;
      if (k < 8) then Notes(note):=Notes(note) + ","; end;
      Notes(note):=Notes(note) + "\n";
    end;
    Notes(note):=Notes(note) + "  }";
    if (j < n) then Notes(note):=Notes(note) + ",\n"; end;
    Notes(note):=Notes(note) + "\n";
  end;
  Notes(note):=Notes(note) + "}"; 
END;

[compsystems]

The icon

[akmon]

Downloaded equation library. I have choosen an easy example:
P=I^2*R
P=V^2/R

Suppose V and R are known variables: V=220 and R is 9 ohm.

This example was able to be solved by MSOLVR fron HP48, because you a have a "concatenated" list of equations, solve P, then solve I.

However the app can´t solve it. Tried various guesses and no way.
Solvesys is able to work as MSOLVR if necesary?

[Han]

(02-14-2017 11:21 PM)akmon Wrote:  Downloaded equation library. I have choosen an easy example:
P=I^2*R
P=V^2/R

Suppose V and R are known variables: V=220 and R is 9 ohm.

This example was able to be solved by MSOLVR fron HP48, because you a have a "concatenated" list of equations, solve P, then solve I.

However the app can´t solve it. Tried various guesses and no way.
Solvesys is able to work as MSOLVR if necesary?

I am not quite sure what you mean by not being able to solve it. If you want it to solve for P and I while fixing the values of V and R, then you will need to tell the solver that the values of V and R are constant. To the far right of the variables are check boxes to set them to be constants (when checked) or variables (when unchecked). If you leave them as variables, then you are actually solving an under-determined system (4 variables with only 2 equations -- there are an infinite number of solutions). On the other hand, setting V and R to be constants gives me the solution P=5377.7778 and I=24.4444 (2 variables, 2 equations -- if a solution exists, it will be unique). I used a guess of 100, 100.

I changed my tolerances because the values supplied (220 and 9) suggest we only need to be within 5 or so decimal places for a reasonable answer (unlike the other problem where the solution itself was in the order of 10^-10 or 10^-11).

Please keep in mind that the conclusion from the app is all based on your tolerances. It stops the iterations if the x_i's are within the tolerance you specify (it's actually not a tolerance, but a check on whether the digits up to a certain order have stopped changing). Then it checks to see if the output values are also within (another) tolerance, concluding a zero if yes. And if not it checks to see if perhaps a minimum (or extremum in general) may have been found. These last two are the only possible, valid conclusions since a "tiny" (or no) change in the x_i's suggests that the gradient is nearly 0 (which occurs at extrema). Here, "tiny" is relative to the tolerance you set. But again, these conclusions are based on your tolerances. So if you set them too large, then it may conclude incorrectly. Set them too small, and you may never see convergence in some situations (not enough precisions to meet the desired tolerance) despite the values actually converging.

That said, I have been considering implementing a hybrid solver just as described by the HP patent. Since the Prime has access to a CAS, it should not be too hard to implement a secondary solver (or maybe this could be the primary solver) in which a formula having only one parameter that is unknown is used to determine the remaining parameter, repeating as much as necessary while cycling through the remaining formulas.

EDIT: Ahh, ok. I downloaded the version listed in the first post, and am seeing the misbehavior. There's a single line in the function ssEvalF that was modified and somehow did not make it into this revision (this is what I get for not using a CVS and developing on multiple copies over multiple machines).

In the meantime, change the line ssCAS(ssCurSys(4,j) + ":=" + vec(j)); to become ssCAS(ssCurSys(4,ssSolveVars(j)) + ":=" + vec(j)); within the function ssEvalF(vec)

I found a few more bugs that I need to iron out in the pinv() routine (for degenerate cases like m x 1 or 1 x n matrices) but I should have an update up soon.

EDIT 2: 0.652 should take care of this, now. SVD was also updated.

[akmon]

Ok, this is another thing. It solves correctly. It´s curious that if I insert R=8, the program founds a zero, but if I input R=9 it says "Bad guess", however the result is correct. The clue is F value. If it is very low you can assume that the solution is quite correct. I´ll make more tests.

By the way, this equation was pending of answer. I´ve read that the app allows to insert complex equations, I tried this: (2+i)*X+2*Y=1+7*i ; (1-i)*X+i*Y=0 but I get a syntax error, don´t know to input the i value.