The following warnings occurred:
Warning [2] count(): Parameter must be an array or an object that implements Countable - Line: 795 - File: showthread.php PHP 7.4.33 (FreeBSD)
File Line Function
/showthread.php 795 errorHandler->error





Post Reply 
Creating an equation library (Updated 03-FEB-2017)
01-09-2015, 09:51 PM (This post was last modified: 03-09-2017 04:25 PM by Han.)
Post: #1
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
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
Creating an equation library (Updated 03-FEB-2017) - Han - 01-09-2015 09:51 PM
RE: Creating an equation library - Han - 01-10-2015, 03:45 AM



User(s) browsing this thread: 5 Guest(s)