Creating an equation library (Updated 03-FEB-2017)
(This post has been edited; the original content of the first post is at the bottom of this post)
MARCH 2017: This project now exists as an actual app. Please see http://hpmuseum.org/forum/thread-7725.html
------------
Version 0.063 now makes use of the SVD package I wrote and posted here: http://www.hpmuseum.org/forum/thread-497...67743.html
This version falls back on the SVD package I wrote for certain classes of matrices -- mostly those with singular values of 0. The Prime's implementation computes the left orthonormal vectors based on the reciprocal of the singular values. So when a singular value is zero, this results in undefined orthonormal vectors. There are still many errors that need to be trapped, but if you stall the SVD program in addition to the code below, you should be able to solve any system, now -- provided your guesses does not result in a point \( \mathbf{x} \) such that \( \mathbf{F}(\mathbf{x} )\) or its Jacobian is undefined.
I will eventually write some documentation, but for now the solve is a simple Newton solver. Briefly,
\[ \mathbf{F}(\mathbf{x}) =
\left[ \begin{array}{c}
f_1 (x_1, x_2, \dotsm, x_n)\\
f_2 (x_1, x_2, \dotsm, x_n) \\
\vdots \\
f_m (x_1, x_2, \dotsm, x_n)
\end{array} \right] \]
where the \( f_i \)'s are the equations set to 0. For example, \( x^2 + xy = 9 \) becomes \( f(x,y) = x^2+xy - 9 \). The linearization of \( \mathbf{F} \)
is
\[ \mathbf{F} (\mathbf{x}_{n+1} ) = \mathbf{F} (\mathbf{x}_n ) + J_\mathbf{F} \cdot \underbrace{(\mathbf{x}_{n+1} - \mathbf{x}_{n} ) }_{\Delta \mathbf{x}}\]
where \( J_{\mathbf{F}} \) is the Jacobian (evaluated at \( \mathbf{x}_{n} \) ). The solver solves for \( \Delta \mathbf{x} \) by computing the pseudo inverse of \( J_\mathbf{F} \) -- which is exactly the inverse if the Jacobian if it is non-singular. I won't get too much SVD and pseudo-inverse until I write up the documentation.
Code:
ssVersion:="Equation Library 0.063 by Han Duong";
// colors
ssFG:=#0h;
ssBG:=#FFFFFFh;
ssR:=#FF0000h;
ssG:=#FF00h;
ssB:=#FFh;
// error messages
ssInvDat:="Invalid library data.";
ssNoDat:="Unable to locate library data.";
ssNullName:="Null name not allowed.";
ssInvName:="Invalid name.";
ssAbort:="Operation aborted; improper initialization.";
ssBadVar:="Conflict with existing variable:\n";
ssSubmitBug:="This should never happen!\nPLease submit bug report.";
ssNotEqn:="Non-equation on page ";
ssNonVar:="No such var. on page ";
ssNoLibDat:="No library data found.\nCreating empty library and empty\ninitial system.";
ssEOW:="Overwrite disabled and name already exists:\n";
ssNotSaved:="SAVE CANCELED!";
ssMustChoose:="Select a system or create a new one.";
ssNoVars:="No variables for this system!";
ssNoSystems:="No systems in library!";
ssBadJac:="Jacobian not invertible; singular or non-square?";
ssFUndef:="Error evaluating F!";
// app messages
ssTNewVar:="Create/Edit Variable";
ssLName:="Name:";
ssLDescr:="Description:";
ssLOWrite:="Overwrite";
ssHName:="Enter the variable name";
ssHDescr:="Enter a description for the variable";
ssHOWrite:="Overwrite if variable already exists?";
ssTSaveSys:="Save Current System?";
ssHSysName:="Enter a name for the system.";
ssLEqns:="Eqns:";
ssLVars:="Vars:";
ssHEqns:="Verify equations in system.";
ssHVars:="Verify variables for system.";
ssTNewSys:="New System";
ssTChooseSys:="Choose a System";
ssHMakeConst:="Make this variable constant?";
ssTDeleteSys:="Select Systems to Delete";
ssHMarkDel:="Mark for deletion?";
ssTSettings:="Equation Library Settings";
ssLSaving:="When switching systems, always";
ssCSaving:={ "prompt to save", "save automatically", "discard all changes"};
ssHSaving:="Select default save behavior";
// misc
ssVLoadEqns:="Populating equations in Symb view...";
ssVLoadVars:="Setting up equation variables...";
ssVCheckVars:="Verifying variables in system...";
ssVCheckLibDat:="Verifying library data...";
ssVUpdateEqns:="Equations modified; updating lib data...";
ssVUpdateVars:="Variables modified; updating lib data...";
ssVNoEqns:="(No equations)";
ssVNoVars:="(No variables)";
ssVVarMade:="Created variable: ";
// View menu options
ssTSEP:="Select Equations Page";
ssTSelPage:="Select New Page";
ssTNewPage:="Create New Page";
ssMsgRow;
export ssInit;
ssLibFN;
export ssCurSys;
export ssCurSysIndex;
export ssCurLib;
ssLibSize;
export ssSysTitles;
ssVarKey:="SSVersion";
export ssEqPages;
export ssEqPage;
export ssNullSys;
export ssWarning;
export ssSolveVars:={};
export ssAutoVar:=1;
// DO NOT MODIFY BELOW HERE
ssGlobal:={ "A", "B", "C", "D", "E", "F", "G",
"H", "I", "J", "K", "L", "M", "N",
"O", "P", "R", "S", "T", "U", "V",
"W", "X", "Y", "Z", "θ",
"Z0", "Z1", "Z2", "Z3", "Z4",
"Z5", "Z6", "Z7", "Z8", "Z9"
};
ssUpdateSys()
begin
local n;
ssEqPages:={1,10};
n:=size(ssCurSys(2));
ssNullSys:=NOT(n);
ssEqPages(1):=ip((n-1)/10);
ssEqPages(2):=(n-1) mod 10;
ssEqPage:=0;
for n from 1 to size(ssCurSys(4)) do
if NOT(ssCurSys(7,n)) then ssSolveVars(0):=ssCurSys(4,n); end;
end;
DelAVars(AVars);
end;
ssInitSSApp()
begin
local j,run:=1;
local cmd:="";
if ssInit then return(0); end;
dimgrob_p(G1,320,240);
ssLibFN:="Equation Library.lib";
ssVarKey:="SSVersion";
ssWarning:=1;
iferr
ssCurLib:=AFiles(ssLibFN);
then
msgbox(ssNoLibDat);
ssCurLib:={};
AFiles(ssLibFN):=ssCurLib;
end;
for j from 0 to 9 do
cmd:="E" + string(j,1,0) + ":=" + string("");
expr(cmd);
end;
// always start app with new system
ssCurSys:={ssTNewSys, {}, {}, {}, {}, {}, {} };
ssLibSize:=size(ssCurLib);
if ssLibSize then
ssSysTitles:=makelist(ssCurLib(j,1),j,1,ssLibSize);
end;
ssCurSysIndex:=0;
ssUpdateSys();
ssInit:=1;
end;
ssPrint(msg)
begin
textout_p(msg,G1,0,12*ssMsgRow,1,ssFG,320,ssBG);
line_p(G1,0,12*ssMsgRow+12,12,12*ssMsgRow+12,ssR);
ssMsgRow:=(ssMsgRow+1) mod 20;
end;
ssError(msg)
begin
textout_p(msg,G1,0,12*ssMsgRow,1,ssR,320,ssBG);
line_p(G1,0,12*ssMsgRow+12,12,12*ssMsgRow+12,ssR);
ssMsgRow:=(ssMsgRow+1) mod 20;
end;
ssWarn(msg)
begin
if ssWarning then
msgbox(msg);
end;
end;
ssCheckType(list,t)
begin
local n:=size(list);
if n then
if type(list(n))<>t then return(1); end;
else
return(1);
end;
return(0);
end;
//******************************************
// Library Data format
// { sys1, sys2, ... , sysN }
//
// Each sys (system) is a list of the form
// {
// "Title",
// { "eq1", "eq2", ... },
// { sel1, sel2, ... },
// { "var1", "var2", ... },
// { "des1", "des2", ... },
// { val1, val2, ... },
// { con1, con2, ... }
// }
//
// eq = equation
// sel = 0 or 1 for selection of corresponding
// equation in the system
// var = variable
// des = description of variable
// val = value of corresponding variable
// con = 0 or 1 for whether var is constant
//
// returns 0 if error or 1 if OK
//******************************************
export ssCheckLibDat()
begin
ssPrint(ssVCheckLibDat);
ssInitSSApp();
local syseqs,j,k;
local lib, isbad:=0;
local types:={2,2,0,2,2,0,0};
iferr
if ssLibSize then
for j from 1 to ssLibSize do
lib:=ssCurLib(j);
if type(lib)<>6 then
ssError(ssInvDat); return(0);
else
k:=size(lib);
if (k<1) OR (k>7) then
ssError(ssInvDat); return(0);
end;
for k from 2 to 7 do
if type(lib(k))<>6 then
ssError(ssInvDat); return(0);
end;
if ssCheckType(lib(k),types(k)) then
ssError(ssInvDat); return(0);
end;
end; // for k
end; // if type(lib)
end; // for j
end; // if ssLibSize
then
ssError(ssNoDat); return(0);
end;
return(1);
end;
// initialize variables
ssInitVars()
begin
local varslist:=ssCurSys(4);
local valslist:=ssCurSys(6);
local n:=size(varslist);
local j,varval;
iferr
if (AVars(ssVarKey) == ssVersion) then
return(0);
end;
then
// nop
end;
ssPrint(ssVLoadVars);
for j from 1 to n do
if pos(ssGlobal,varslist(j)) then
expr(varslist(j) + ":=" + valslist(j));
else
varval:=CAS(EVAL(varslist(j)));
purge(EVAL(varslist(j)));
if (type(CAS(EVAL(varslist(j)))) <> 8) then
ssWarn(ssBadVar+varslist(j));
end;
if ( (type(varval) == 0) OR (type(varval) == 3) ) then
AVars(varslist(j)):=varval;
else
AVars(varslist(j)):=valslist(j);
end;
end;
end;
AVars(ssVarKey):=ssVersion;
end;
// checks to see if AVars and system of equations
// vars match; excludes version key
export ssCheckVars()
begin
ssPrint(ssVCheckVars);
ssInitSSApp();
ssInitVars();
local j,n;
local avarlist:=AVars();
local varlist:=ssCurSys(4);
n:=size(varlist);
for j from 1 to n do
if NOT(pos(avarlist,varlist(j)) OR pos(ssGlobal,varlist(j))) then
return(0);
end;
end;
n:=size(avarlist);
for j from 1 to n do
if ( (avarlist(j) <> ssVarKey) AND NOT(pos(varlist,avarlist(j))) ) then
return(0);
end;
end;
return(1);
end;
// do we have any equations
// 1: yes, equations in symb view
// 0: no equations at all
// -1: yes, but only saved in library
ssHaveEqns()
begin
local j;
for j from 0 to 9 do
if ISCHECK(j) then return(1); end;
end;
if (ssLibSize == 0) then return(0); end;
if (type(ssCurSys) <> 6) then return(0); end;
if size(ssCurSys) then
if size(ssCurSys(2)) then
return(-1);
else
return(0);
end;
else
return(0);
end;
end;
// UI for creating new variables
ssCreateVar()
begin
local run:=1;
local n;
local varname:="", varinfo:="";
local ow:=0;
ssInitSSApp();
ssInitVars();
while run do
if input(
{
{varname, [2], {30, 65, 0}},
{varinfo, [2], {30, 65, 1}},
{ow, 1, {30,10,2}}
},
ssTNewVar,
{ ssLName, ssLDescr, ssLOWrite },
{ ssHName, ssHDescr, ssHOWrite }
)
then
if size(varname)<1 then
msgbox(ssNullName);
else
n:=pos(AVars,varname);
if n then
if ow then
n:=pos(ssCurSys(4),varname);
ssCurSys(5,n):=varinfo;
run:=0;
else
msgbox(ssEOW + varname);
end;
else
iferr
if NOT(pos(ssGlobal,varname)) then
AVars(varname):=0;
end;
then
msgbox(ssInvName);
else
purge(EVAL(varname));
ssCurSys(4,0):=varname;
ssCurSys(5,0):=varinfo;
ssCurSys(6,0):=0;
ssCurSys(7,0):=0;
msgbox(ssVVarMade + varname);
varname:=""; varinfo:="";
end;
end;
end; // end null name
else
run:=0;
end;
end; // end while
end;
// save current equations
ssSaveEqns()
begin
local j,k,n1,n2;
local cmd:="";
local oldeqn:="";
local neweqn:="";
local markadd:=0;
local neweqns:={};
local teqnlist:={};
local tsellist:={};
if ssCheckVars() then
ssPrint(ssVLoadEqns);
n1:=ssEqPage*10+1;
if (ssEqPage < ssEqPages(1)) then
n2:=n1+9;
else
n2:=n1+ssEqPages(2);
end;
for j from n1 to n2 do
k:= j mod 10;
cmd:="string(E" + string(k,1,0) + ",1,12)";
iferr
neweqn:=expr(cmd);
then
neweqn:="";
else
if ssNullSys then
oldeqn:="";
else
cmd:="E" + string(k,1,0) + ":=" + string(ssCurSys(2,j));
iferr expr(cmd); then end;
cmd:="string(E" + string(k,1,0) + ",1,12)";
iferr oldeqn:=expr(cmd); then end;
end;
if (oldeqn <> neweqn) then
if ssNullSys then
ssCurSys(2,0):=neweqn;
else
ssCurSys(2,j):=neweqn;
end;
end;
end;
cmd:="E" + string(k,1,0) + ":=" + string(neweqn);
iferr expr(cmd); then
// msgbox(ssSubmitBug);
end;
if ISCHECK(k) then
ssCurSys(3,j):=1;
else
ssCurSys(3,j):=0;
end;
if (size(neweqn) == 0) then
ssCurSys(3,j):=-1; // mark for delete
end;
end; // end for j
n2:=n1+9;
for k from j to n2 do
n1:=k mod 10;
cmd:="string(E" + string(n1,1,0) + ",1,12)";
iferr
neweqn:=expr(cmd);
then
// nop
else
if size(neweqn) then
neweqns(0):=neweqn;
neweqns(0):=ISCHECK(n1);
markadd:=1;
end;
end;
end; // end for k
// rebuilt existing equations
k:=0;
n2:=size(ssCurSys(2));
for j from 1 to n2 do
if (ssCurSys(3,j) >= 0) then
teqnlist(0):=ssCurSys(2,j);
tsellist(0):=ssCurSys(3,j);
k:=k+1;
end;
end;
ssCurSys(2):=teqnlist;
ssCurSys(3):=tsellist;
if markadd then
n2:=size(neweqns)/2;
for j from 1 to n2 do
ssCurSys(2,0):=neweqns(2*j-1);
ssCurSys(3,0):=neweqns(2*j);
k:=k+1;
end;
end;
ssEqPages(1):=ip((k-1)/10);
ssEqPages(2):=(k-1) mod 10;
ssEqPage:=min(max(0,ssEqPage),ssEqPages(1));
if k then
ssNullSys:=0;
else
ssNullSys:=1;
end;
return(1);
else
msgbox(ssAbort);
return(0);
end;
end;
// set up equations in Symb view
export ssLoadSymb()
begin
local j,k,n1,n2;
local cmd:="";
// no equations to populate Symb view?
if ssNullSys then
for j from 0 to 9 do
cmd:="E" + string(j,1,0) + ":=" + string("");
expr(cmd);
end;
return(1);
end;
ssPrint(ssVLoadEqns);
n1:=ssEqPage*10+1;
if (ssEqPage < ssEqPages(1)) then
n2:=n1+9;
else
n2:=n1+ssEqPages(2);
end;
for j from n1 to n2 do
k:= j mod 10;
cmd:="E" + string(k,1,0) + ":=" + string(ssCurSys(2,j));
iferr expr(cmd); then end;
if (ssCurSys(3,j)==1) then CHECK(k); else UNCHECK(k); end;
end;
n2:=n1+9;
for k from j to n2 do
n1:=k mod 10;
cmd:="E" + string(n1,1,0) + ":=" + string("");
expr(cmd);
end;
return(1);
end;
// select an equations page
ssSelectEqPage()
begin
ssInitSSApp();
local j;
local n:=ssEqPages(1)+1;
local page:=ssEqPage+1;
local curpage:=page;
local pagelist:={};
local title:="[" + string(page,1,0) + "/" + string(ssEqPages(1)+1,1,0) + "] " + ssTSelPage;
pagelist:=makelist(j,j,1,n);
pagelist(0):=ssTNewPage;
if choose(page,title,pagelist) then
if (curpage <> page) then
ssSaveEqns();
if (page == (n+1)) then
if (ssNullSys == 0) then
ssEqPages(1):=n+1;
// fill last page
if (ssEqPages(2) < 9) then
n:=9-ssEqPages(2);
for j from 1 to n do
ssCurSys(2,0):="";
ssCurSys(3,0):=-1; // mark for deletion
end;
end;
end;
// create new page
for j from 1 to 10 do
ssCurSys(2,0):="";
ssCurSys(3,0):=-1; // mark for deletion
end;
n:=size(ssCurSys(2));
ssEqPage:=ip((n-1)/10);
ssEqPages(1):=ssEqPage;
ssEqPages(2):=(n-1) mod 10;
ssNullSys:=0;
else
// not creating new blank page
ssEqPage:=min(max(0,page-1),ssEqPages(1));
end; // if newpage or existing page
end; // if selpage <> page
ssLoadSymb();
end; // if choose
end;
// set new system
export ssSetSystem(n)
begin
if (n > ssLibSize) then
return(0);
end;
ssCurSysIndex:=n;
ssCurSys:=ssCurLib(n);
ssUpdateSys();
ssInitVars();
ssLoadSymb();
end;
ssErrorPage(n,msg,eqn)
begin
local p:=ip((n-1)/10)+1;
msgbox(msg + n + ":\n" + eqn);
end;
// automatically add existing variable to vars list
// should only happen with global variables
ssAutoAddVar(var)
begin
local val;
local t;
iferr
val:=expr(var);
then
val:=0;
end;
t:=type(val);
if ( (t<>0) OR (t<>3) ) then
val:=0;
end;
// convert existing vars to app var
// exceptions: global system vars cannot
// be app vars in current firmware
if NOT(pos(ssGlobal,var)) then
AVars(var):=val;
end;
ssCurSys(4,0):=var;
ssCurSys(5,0):="";
ssCurSys(6,0):=val;
ssCurSys(7,0):=0;
end;
// check if equations are valid
// assumes equations have been saved
export ssCheckEqns()
begin
local n:=size(ssCurSys(2));
local eqn:="";
local j,k,v;
local varlist:={};
local var:="";
if n then
// check each equation
for j from 1 to n do
eqn:=ssCurSys(2,j);
// do we have = symbol?
if instring(eqn,"=") then
iferr eqn:=string(E0,1,12); then eqn:=""; end;
E0:=ssCurSys(2,j);
varlist:=LNAME(E0);
v:=size(varlist);
if v then
for k from 1 to v do
var:=string(varlist(k));
if (pos(ssCurSys(4),var) == 0) then
if ssAutoVar then
ssAutoAddVar(var);
else
ssErrorPage(ip((j-1)/10)+1,ssNonVar,ssCurSys(2,j)+"\n> "+varlist(k));
iferr E0:=eqn; then end;
return(0);
end;
end;
end;
else
// constants only? not equation
ssErrorPage(j,ssNotEqn,ssCurSys(2,j));
iferr E0:=eqn; then end;
return(0);
end;
iferr E0:=eqn; then end;
else
// no equal symbol; report page
ssErrorPage(j,ssNotEqn,eqn);
return(0);
end;
end;
// all eqns seem ok
return(1);
else
// no equations to check
return(-1);
end;
end;
// save current system
// no null names but does not check
// validity of names beyond whether
// AVars will accept such a name
ssSaveSystem()
begin
ssSaveEqns();
local sysname:=ssCurSys(1);
local eqns:=ssCurSys(2);
local vars:=ssCurSys(4);
local eqn,var,n;
local run:=1;
// bug in input when choose lists are empty
if (size(eqns) < 1) then eqns:={ ssVNoEqns }; end;
if (size(vars) < 1) then vars:={ ssVNoVars }; end;
while run do
if input(
{
{ sysname, [2], { 15, 75, 0 } },
{ eqn, eqns, { 15, 75, 1 } },
{ var, vars, { 15, 75, 2 } }
},
ssTSaveSys,
{ ssLName, ssLEqns, ssLVars },
{ ssHSysName, ssHEqns, ssHVars }
)
then
if size(sysname) then
// maybe later add check for dupes
n:=pos(ssSysTitles, sysname);
if ssCheckEqns() then
ssCurSys(1):=sysname;
ssCurLib(ssCurSysIndex):=ssCurSys;
AFiles(ssLibFN):=ssCurLib;
if (ssCurSysIndex == 0) then ssCurSysIndex:=size(ssCurLib); end;
ssLibSize:=size(ssCurLib);
else
ssWarn(ssNotSaved);
end;
run:=0;
else
msgbox(ssNullName);
end;
else
// user canceled
ssWarn(ssNotSaved);
run:=0;
end;
end; // end while;
end;
// create a system
ssNewSystem()
begin
ssCurSys:={ssTNewSys, {}, {}, {}, {}, {}, {} };
ssCurSysIndex:=0;
ssUpdateSys();
ssInitVars();
ssLoadSymb();
startview(0,1);
end;
ssChooseSys()
begin
local n,t;
local run:=1;
local systitles:={};
n:=ssCurSysIndex;
iferr
ssCurLib:=AFiles(ssLibFN);
then
ssWarn(ssNoLibDat);
ssCurLib:={};
AFiles(ssLibFN):=ssCurLib;
end;
ssLibSize:=size(ssCurLib);
if ssLibSize then
systitles:=makelist(ssCurLib(n,1),n,1,ssLibSize);
else
systitles:={};
end;
ssSysTitles:=systitles;
systitles(0):=ssTNewSys;
t:=choose(n, ssTChooseSys, systitles);
if t then
if ( n <= ssLibSize ) then
ssSetSystem(n);
else
ssNewSystem();
end;
return(1);
end;
return(0);
end;
// UI for deleting systems
export ssDeleteSys()
begin
local cmd:="input({";
local systitles;
local j;
local selected:={};
local lib:={};
if ssLibSize then
selected:=makelist(0,j,1,ssLibSize);
systitles:=makelist(ssCurLib(j,1),j,1,ssLibSize);
for j from 1 to ssLibSize do
cmd:=cmd + "{selected(" + string(j,1,0) + "),0,{94,5," + string(j-1,1,0) + "}}";
if (j<ssLibSize) then cmd:=cmd + ",\n"; end;
end;
cmd:=cmd + "},\n" + string(ssTDeleteSys) + ",\n{";
for j from 1 to ssLibSize do
cmd:=cmd + string(ssCurLib(j,1));
if (j<ssLibSize) then cmd:=cmd + ",\n"; end;
end;
cmd:=cmd + "},\n{";
for j from 1 to ssLibSize do
cmd:=cmd + string(ssHMarkDel);
if (j<ssLibSize) then cmd:=cmd + ",\n"; end;
end;
cmd:=cmd + "},\n" + string(selected,1,0) + ",\n" + string(selected,1,0) + ")";
Notes("eqlib_debug.log"):=cmd;
j:=expr(cmd);
if j then
// current system deleted? if so set index to 0
if ssCurSysIndex then
if selected(ssCurSysIndex) then
ssCurSysIndex:=0;
end;
end;
// rebuild library and adjust index
for j from 1 to ssLibSize do
if NOT(selected(j)) then
lib(0):=ssCurLib(j);
if (j < ssCurSysIndex) then
ssCurSysIndex:= ssCurSysIndex - 1;
end;
end;
end;
ssCurLib:=lib;
AFiles(ssLibFN):=ssCurLib;
end;
else
ssWarn(ssNoSystems);
end;
end;
// variable browser
export ssVarBrowser()
begin
local vals:=ssCurSys(6);
local cons:=ssCurSys(7);
local cmd:="input({";
local n:=size(vals);
local j,k;
local tmp:="";
local xloc:="15,75";
if (n<1) then
msgbox(ssNoVars);
return(0);
end;
k:=1;
for j from 1 to n do
k:=max(k,size(ssCurSys(4,j)));
end;
if (k > 3) then xloc:="30,60"; end;
// build input( ... )
// { variable, type, { pos } }
for j from 1 to n do
k:=pos(ssGlobal, ssCurSys(4,j));
if (k<28) then tmp:="0"; else tmp:="0,3"; end;
cmd:=cmd + "{vals(" + string(j,1,0) + "),[" + tmp + "],{" + xloc + "," + string(j-1,1,0) + "}}, ";
cmd:=cmd + "{cons(" + string(j,1,0) + "),0,{94,5," + string(j-1,1,0) + "}}";
if ((n>1) AND (j<n)) then cmd:=cmd + ",\n"; end;
end;
// add title
cmd:=cmd + "},\n" + string(ssCurSys(1) + " Variables") + ",\n{";
// add labels
for j from 1 to n do
cmd:=cmd + string(ssCurSys(4,j) + "=") + ", " + string("");
if ((n>1) AND (j<n)) then cmd:=cmd + ",\n"; end;
end;
cmd:=cmd + "},\n{";
// add help
for j from 1 to n do
cmd:=cmd + string(ssCurSys(4,j) + ": " + ssCurSys(5,j)) + "," + string(ssHMakeConst);
if ((n>1) AND (j<n)) then cmd:=cmd + ",\n"; end;
end;
cmd:=cmd + "},\n";
// add init and reset vals
tmp:="{";
for j from 1 to n do
tmp:=tmp + vals(j) + "," + cons(j);
if ((n>1) AND (j<n)) then tmp:=tmp + ","; end;
end;
tmp:=tmp + "}";
cmd:=cmd + tmp + ",\n" + tmp + ")";
Notes("eqlib_debug.log"):=cmd;
j:=expr(cmd);
if j then
ssCurSys(6):=vals;
ssCurSys(7):=cons;
for j from 1 to n do
if pos(ssGlobal,ssCurSys(4,j)) then
expr(ssCurSys(4,j) + ":=" + vals(j));
else
AVars(ssCurSys(4,j)):=vals(j);
end;
end;
end;
return(j);
end;
// build F, F(var1, var2, ...) and Jacobian
// assumes all equations and variables have been
// validated; also assumes equations are saved and
// that selectd equations form a non-null list
// saves F as ssF (CAS variable)
export casMakeF()
begin
local cmd:="['";
local cmd2:="[";
local j,n;
local varlist:={};
local eqnlist:={};
local k:=size(ssCurSys(4));
local vstr1:="";
local vstr2:="";
n:=size(ssCurSys(2));
for j from 1 to n do
if ssCurSys(3,j) then
eqnlist(0):=ssCurSys(2,j);
end;
end;
for j from 1 to k do
if NOT(ssCurSys(7,j)) then
varlist(0):=j;
end;
end;
ssSolveVars:=varlist;
n:=k-1;
for j from 1 to n do
vstr1:=vstr1 + ssCurSys(4,j) + ",";
end;
vstr1:=vstr1 + ssCurSys(4,k);
n:=size(varlist)-1;
for j from 1 to n do
vstr2:=vstr2 + ssCurSys(4,varlist(j)) + ",";
end;
vstr2:=vstr2 + ssCurSys(4,varlist(n+1));
// ssF:=( vars ) -> [ f1, f2, ... ]
n:=size(eqnlist)-1;
for j from 1 to n do
cmd:=cmd + eqnlist(j) + ")','";
cmd2:=cmd2 + eqnlist(j) + "),";
end;
cmd:=cmd + eqnlist(n+1) + ")']";
cmd:=replace(cmd,"=","-(");
CAS(EVAL("ssF:=" + cmd));
cmd2:=cmd2 + eqnlist(n+1) + ")]";
cmd2:=replace(cmd2,"=","-(");
CAS(EVAL("ssFx:=(" + vstr1 + ")->(" + cmd2 + ")"));
cmd:="ssJ:=transpose(diff(" + cmd + ",[" + vstr2 + "]))";
CAS(EVAL(cmd));
end;
// pseudo inverse using SVD
export ssPInv()
begin
local l:=CAS("svd(ssJ)");
local r:=size(l(2)); r:=r(1);
local j,n,ns,s;
local tol:=1E-30;
local d:=[[0]];
// built-in SVD fails on some singular matrices
// failed SVD matrices returned as list of lists
if ((type(l(1)) == 6) OR (type(l(3)) == 6)) then
d:=CAS("pinv(ssJ,1e-30)");
return(d);
end;
// sort sing. vals from largest to smallest
// and adjust col/row accordingly
for n from r downto 2 do
ns:=1;
s:=l(2,1);
for j from 2 to n do
if (l(2,j) <= s) then
ns:=j; s:=l(2,j);
end;
end;
if (ns < n) then
l(2,ns):=l(2,n);
l(2,n):=s;
d:=l(1);
l(1):=swapcol(d,ns,n);
d:=l(3);
l(3):=swapcol(d,ns,n);
end;
end;
for n from 1 to r do
s:=l(2,n);
if (s >= tol) then
l(2,n):=1/s;
else
l(2,n):=0;
end;
end;
j:=size(l(3)); n:=j(1);
j:=size(l(1)); r:=j(1);
d:=makemat(ifte(I==J,l(2,I),0),n,r);
// make sure to use CAS.trn for Hermetian
return(l(3)*d*CAS.trn(EVAL(l(1))));
end;
// Newton's Method
// dx = - J^(-1) * F
// x[n+1] = x[n] + dx
export casNewton()
begin
local j,k,nsv,nv;
local var;
local tolx:=1.0E-8;
local tolF:=1.0E-10; // exit if F.F < tolF
local nsF=2*tolF; // initial squared norm of F
local vecNewt:={};
// L1:={}; // for debugging only
local pinv:=[[0]];
local dx:=2*tolx;
nv:=size(ssCurSys(4));
nsv:=size(ssSolveVars);
for j from 1 to nsv do
vecNewt(0):=ssCurSys(6,ssSolveVars(j));
// L1(0):=ssCurSys(6,ssSolveVars(j));
end;
for k from 1 to 100 do
nsF:=CAS("approx(sqrt(ssF*ssF))");
if (type(nsF) == 2) then
msgbox(ssFUndef);
break;
end;
vecNewt(0):=nsF;
if (k == 1) then
vecNewt(0):=0;
else
vecNewt(0):=abs(dx);
end;
if (nsF < tolF) then break; end;
if (abs(dx) < tolx) then break; end;
iferr
pinv:=ssPInv();
dx:=-pinv*CAS.transpose("approx(ssF)");
then
msgbox(ssBadJac);
break;
end;
for j from 1 to nsv do
var:=ssCurSys(4,ssSolveVars(j));
if pos(ssGlobal,var) then
expr(var + ":=" + var + "+" + dx(j,1));
else
AVars(var):=AVars(var) + dx(j,1);
end;
ssCurSys(6,ssSolveVars(j)):=expr(var);
end; // for j
for j from 1 to nsv do
vecNewt(0):=ssCurSys(6,ssSolveVars(j));
// L1(0):=ssCurSys(6,ssSolveVars(j));
end;
end; // for k
return(CAS.list2mat(EVAL(vecNewt),EVAL(nsv+2)));
end;
START()
begin
// startview(-1,1);
ssInitSSApp();
ssInitVars();
startview(0,1);
end;
// Equations Editor
Symb()
begin
if ssInit then
ssSaveEqns();
else
ssInitSSApp();
ssInitVars();
end;
ssLoadSymb();
startview(0,1);
end;
// System of Equations browswer
Plot()
begin
if ssInit then
ssSaveSystem();
else
ssInitSSApp();
ssInitVars();
end;
ssChooseSys();
startview(0,1);
end;
// run solver
Num()
begin
local matNewt:=[[0]];
local j;
local tlist:={};
local labels:={};
if ssInit then
ssSaveEqns();
if NOT(ssCheckEqns()) then
return(0);
end;
if NOT(ssHaveEqns()) then
return(0); // exit
end;
if ssVarBrowser() then
casMakeF();
matNewt:=casNewton();
tlist(0):=ssCurSys(1);
tlist(0):={};
for j from 1 to size(ssSolveVars) do
labels(0):=ssCurSys(4,ssSolveVars(j));
end;
labels(0):="ǁFǁ";
labels(0):="ǁΔxǁ";
tlist(0):=labels;
editmat(matNewt,tlist);
end;
else
ssInitSSApp();
ssInitVars();
end;
end;
// Set up initial guess
NumSetup()
begin
if ssInit then ssSaveEqns(); end;
ssVarBrowser();
end;
// select equations page (View Menu)
view "Select Page", ssMenuEqnPage()
begin
ssSelectEqPage();
startview(0,1);
end;
// UI for creating new variables (View Menu)
view "Add/Edit Variable", ssMenuNewVar()
begin
ssCreateVar();
end;
// reload saved equations (View Menu)
view "Restore Equations", ssMenuLoadEqns()
begin
if ssCurSysIndex then
ssCurSys:=ssCurLib(ssCurSysIndex);
ssUpdateSys();
ssInitVars();
end;
ssLoadSymb();
startview(0,1);
end;
// save current system (View Menu)
view "Save System", ssMenuSaveSys()
begin
ssSaveSystem();
end;
// create new system (View Menu)
view "New System", ssMenuNewSys()
begin
ssSaveSystem();
ssNewSystem();
end;
// delete existing system (View Menu)
view "Delete System", ssMenuDelSys()
begin
ssSaveSystem();
ssDeleteSys();
end;
// select system (View Menu)
view "Select System", ssMenuSetSys()
begin
ssSaveSystem();
if ssChooseSys() then startview(0,1); end;
end;
Original first post below -- info is out of date due to revisions
-------------------------------------------------------------------------------------------
Here's a brief outline of my ideas.
1) Each category of equations will be a list of lists of the following form:
Code:
{
"Title",
{ "equation1", "equation2", ..., "equationn" },
{ "var1", "var2", ..., "varn" },
{ "vardescription1", "vardescription2", ..., "vardescriptionn" }
}
Perhaps later on we can include graphics by means of lists of point definitions, line definitions, and polygon definitions. Ideally, the user would not have to specify what the variable list is, and the app will be smart enough to parse through each equation string to identify what the variables are.
2) A choose screen will allow users to select which category of equations to work with
3) Once a category is selected, all the equations will be shown in an INPUT() window with check boxes next to them. The user then checks the equations to use.
4) The app then creates another INPUT() window for variables and their initial values. Each variable may also be checked to specify whether the quantity should be treated as a variable that may change during the solve process, or if it is a constant (much like in the SOLVESYS library for the HP48 and HP49).
5) The equations are then modified to insert constants where appropriate, and then call fsolve() to solve the system.
This is the general scheme. What would also be nice is an interface for users to add in new categories of equations. Anyway, this is something for a later update. For now, I'll start working on building the initial interface with the assumption that categories of equations already exist. This is where, if you want to contribute, you can help by posting the categories as shown in the format above.
I'll start with a simple example:
Code:
{
"Linear Motion",
{
"x=x0+v0*t+1/2*a*t^2",
"x=x0+v*t-1/2*a*t^2",
"x=x0+1/2*(v0+v)*t",
"v=v0+a*t"
},
{ "x", "x0", "v0", "v", "t", "a" },
{
"final position",
"initial position",
"initial velocity",
"final velocity",
"time",
"acceleration"
}
}
Graph 3D | QPI | SolveSys
|
Search
Member List
Calendar
Help
/
Creating an equation library (Updated 03-FEB-2017)
