Creating an equation library (Updated 03-FEB-2017)

[Han]

So the equation library was put on hold until I could get a working SVD that was a little more robust than the built-in one. The one suggested by compsystems is not numerically stable. We use the SVD to compute the pseudo-inverse of the Jacobian matrix in Newton's method, and then a backtracking linesearch to find a suitable step size. Solving \( J \cdot \Delta \mathbf{x} = - \mathbf{F} \) can be done (quicker) using other matrix factorizations (such as QR or LU) and perhaps I will add the option to choose later on. For now, though, the pseudo-inverse from the SVD seems the most robust, albeit the slowest of them all.

There is also a backtracking linesearch as implemented in SOLVESYS from the HP48 days. One can read more about it from Numerical Recipes -- The Art of Computer Programming, 3rd ed, sections 9.6 and 9.7. The messages from root finding have the same meaning as in SOLVESYS (e.g. "Zero" means the solver found an \( \mathbf{x} \) such that \( \lVert \mathbf{F} \rVert_2 < \epsilon \) for some suitable tolerance \( \epsilon \). These messages need further fine tuning in cases of ill conditioned matrices and ill posed problems.

I do not intend to keep the current interface -- it's too clunky and just does not look as nice as it could. There are also plans to switch over to app format so as to take advantage of AFiles for data storage (e.g. storing the equations and possibly pictures/diagrams) so that everything is kept together. Right now, all equations must be saved in a Note name "EqLib" (see earlier posts).

Here's the current work in progress for those who want to tinker with an equation library / solvesys program. In order for this to work, please read earlier posts on setting up the "library data" in a note named "EqLib" and make sure to also install the SVD program found here: http://www.hpmuseum.org/forum/thread-4976.html

Code:

#pragma mode( separator(.,;) integer(h32) )
eql_MAXITS:=200;
eql_TOLF:=1.0e-5;
eql_STPMX:=100;
eql_TOLX:=1.0e-4;
eql_TOLLSQ:=1.0e-4;
eql_SVDEPS:=1.0e-10;

eql_data;
eql_varlist;
eql_eqnslist;
eql_initvals;
eql_results;

vecF:=[0];
oldvecF:=[0];
vecX:=[0];
oldvecX:=[0];
vecD:=[0];
n1, n2;
a1, a2;
maxstep,delta;
numeq,numvars;
tjc;
costol;
vecgrad;

// this is for debugging
export inputstr;


eql_init();
eql_cmenu();
eql_setup();
eql_dovars();
eql_dosolve();
eql_newt();
eql_iterate();
eql_checktol();

export eqlib()
begin
  eql_init();
  local c:=eql_cmenu();
  local tmp;

  if c then
    tmp:=eql_setup(c);
  end;

  if tmp then
    tmp:=eql_dovars(c,1);
  end;

  if tmp then
    repeat
      eql_dosolve(c);
      tmp:=eql_dovars(c,0);
    until tmp==0;
  end;

end;

eql_init()
begin
  local tmp;

  iferr tmp:=Notes("EqLib") then
    msgbox("No equation library data!");
  else
    eql_data:=expr(tmp);
  end;
end;

eql_cmenu()
begin
  local c;
  local n:=size(eql_data);
  local categories:=makelist(eql_data(X,1),X,1,n);

  choose(c,"Equation Library", categories);

  return(c);
end;

eql_setup(c)
begin
  local eqns:=eql_data(c);
  local tmp:="input({";
  local n:=size(eqns(2));
  local i;
  eql_eqnslist:=makelist(1,X,1,n);

  for i from 1 to n do
    tmp:=tmp+"{eql_eqnslist(" + string(i,1,12) + "),0,{94,5," + string(i-1,1,12) + "}}";
    if (n>1) AND (i<n) then tmp:=tmp + ","; end;
  end;
  tmp:=tmp + "}," + string(eqns(1)) + ",{";

  for i from 1 to n do
    tmp:=tmp + string(eqns(2,i));
    if (n>1) AND (i<n) then tmp:=tmp + ","; end;
  end;
  tmp:=tmp + "},{";

  for i from 1 to n do
    tmp:=tmp + string("Use this equation?");
    if (n>1) AND (i<n) then tmp:=tmp + ","; end;
  end;
  tmp:=tmp + "})";
  n:=expr(tmp);
  return(n);
end;


// determine if variable currently
// has a value as viewable from CAS
// and configure constants
eql_dovars(c,f)
begin
  local eqns:=eql_data(c);
  local tmp:="";
  local n:=size(eqns(3));
  local i; local t;
  if f then
    eql_varlist:=makelist(0,X,1,n);
    eql_initvals:=makelist(0,X,1,n);
    for i from 1 to n do
      tmp:=string(CAS(eqns(3,i)));
      if tmp<>eqns(3,i) then
        eql_initvals(i):=CAS(eqns(3,i));
      else
        eql_initvals(i):=0;
      end;
    end;
  end;

  tmp:="input({";

  for i from 1 to n do
    tmp:=tmp + "{eql_initvals(" + string(i,1,12) + "),[0],{30,50," + string(i-1,1,12) + "}},";
    tmp:=tmp + "{eql_varlist(" + string(i,1,12) + "),0,{94,5," + string(i-1,1,12) + "}}";
    if (n>1) AND (i<n) then tmp:=tmp + ","; end;
  end;
  tmp:=tmp + "}," + string(eqns(1) + " Variables") + ",{";

  for i from 1 to n do
    tmp:=tmp + string(eqns(3,i) + "=") + "," + string("");
    if (n>1) AND (i<n) then tmp:=tmp + ","; end;
  end;
  tmp:=tmp + "},{";

  for i from 1 to n do
    tmp:=tmp + string("Enter the value for the " + eqns(4,i)) + ",";
    tmp:=tmp + string("Make this variable a constant?");
    if (n>1) AND (i<n) then tmp:=tmp + ","; end;
  end;
  tmp:=tmp + "},";

  local tmp2:="{";

  for i from 1 to n do
    tmp2:=tmp2 + eql_initvals(i) + "," + eql_varlist(i);
    if (n>1) AND (i<n) then tmp2:=tmp2 + ","; end;
  end;
  tmp2:=tmp2 + "}";

  tmp:=tmp + tmp2 + "," + tmp2 + ")";
  t:=expr(tmp);

  if t then
    for i from 1 to n do
      if eql_varlist(i) then
        tmp:=eqns(3,i) + ":=" + eql_initvals(i);
        CAS(tmp);
      else
        purge(eqns(3,i));
      end;
    end;
  end;
  
  return(t);
end;

eql_dosolve(c)
begin
  local eqns:=eql_data(c);
  local i,j,k,s;
  local tmp:="eql_tjacobian:=diff(";
  local tjc;
  local eqnstr:="[";
  local varstr:="[";

  numeq:=ΣLIST(eql_eqnslist);
  numvars:=size(eqns(3))-ΣLIST(eql_varlist);
  vecX:=[0];
  vecX(numvars):=0; oldvecX:=vecX; vecD:=vecX;
  vecF:=[0];
  vecF(numeq):=0; oldvecF:=vecF;

  // create "[ eq1, eq2, ..., eqn ]"
  k:=numeq; s:=size(eqns(2));
  for i from 1 to s do
    if eql_eqnslist(i) then
      k:=k-1;
      eqnstr:=eqnstr + replace(eqns(2,i),"=","-(") + ")";
      if k then eqnstr:=eqnstr + ","; end;
    end;
  end;
  eqnstr:=eqnstr + "]";

  eql_results:=CAS(eqnstr);
  s:=size(eql_results);
  for i from 1 to s do
    if type(eql_results(i))==0 then msgbox("All constants?"); return(0); end;
  end;

  // create "[ var1, var2, ..., varn ]" and set up initial guess
  k:=numvars;
  if k==0 then msgbox("No variables to solve for!"); return(0); end;
  vecX(k):=0; j:=1; s:=size(eqns(3));
  for i from 1 to s do
    if eql_varlist(i)==0 then
      k:=k-1; vecX(j):=eql_initvals(i); j:=j+1;
      varstr:=varstr + eqns(3,i);
      if k then varstr:=varstr + ","; end;
    end;
  end;
  varstr:=varstr + "]";

  tmp:=tmp + eqnstr + "," + varstr +")";
  CAS(tmp);

  maxstep:=max(abs(vecX),numvars)*eql_STPMX;
  for i from 1 to eql_MAXITS do
    eql_newt(eqnstr,varstr);
    tmp:=eql_checktol();

    if tmp(2) then
      tmp:=varstr + ":=" + vecX;
      CAS(tmp);

      k:=1; s:=size(eql_varlist);
      for j from 1 to s do
        if eql_varlist(j)==0 then
          eql_initvals(j):=vecX(k); k:=k+1;
        end;
      end;

      return;
    end;

  end;
  msgbox("Max iterations exceeded; bad guess?");

end;

eql_newt(eqnstr,varstr)
begin
  local m1,m2,t1,t2;
  local a,b,d,tmp;
  local jac,eps;

  tmp:=varstr + "=" + vecX;
//  tmp:=replace(tmp, "ᴇ", "E");
  tjc:=CAS.subst("eql_tjacobian",tmp);
  vecF:=CAS.subst(eqnstr,tmp);
  vecgrad:=tjc*vecF;
  tjc:=trn(tjc);

  // hack for CAS/Home issues
  jac:=tjc; eps:=eql_SVDEPS;

  // use pseudo inv to solve J*(xnew-xold)=-F
  vecD:=pinv(jac,eps)*vecF;
  oldvecF:=vecF; oldvecX:=vecX;
  n2:=.5*dot(vecF,vecF);

  tmp:=abs(vecD);
  if tmp>maxstep then vecD:=maxstep/tmp * vecD; end;

  a1:=1;

  // linesearch
  while eql_iterate(eqnstr,varstr) do
    if a1==1 then
      a:=1/(n1/n2 + 1);
    else
      t1:=n1+n2*(2*a1-1);
      t2:=m1-m2+2*a2*n2;
      a:=1/(a1-a2)*(t1/a1^2 - t2/a2^2);
      b:=-a2*t1/a1^2 + a1*t2/a2^2;
      if a==0 then
        a:=n2/b;
      else
        d:=b*b+6*a*n2;
        if d<0 then
          a:=0.5*a1;
        else
          if b<=0 then
            a:=(-b + sqrt(d)) / (3*a);
          else
            a:=-2*n2/(b+sqrt(d));
          end;
        end;
      end;
    end;
    a2:=a1; m1:=n1; m2:=n2;
    a1:=max(.1*a1, min(a,.5*a1));
  end;

end;

eql_iterate(eqnstr,varstr)
begin
  local tmp;

  vecX:=oldvecX-a1*vecD;
  delta:=abs(a1*vecD)/max(1,abs(vecX));
  tmp:=varstr + "=" + vecX;
//  tmp:=replace(tmp, "ᴇ", "E");
  vecF:=CAS.subst(eqnstr,tmp);
  n1:=0.5*dot(vecF,vecF);

  if n1<eql_TOLF then return(0); end;
  if n1==n2 then return(0); end;
  if a1<1E-10 then return(0); end;
  return(n1>n2);

end;

eql_checktol()
begin
  if n1<eql_TOLF then msgbox("Zero"); return({1,1}); end;
  if (n1==n2) AND (delta > eql_TOLX) then msgbox("Bad guess(es)"); return({1,0}); end;
  if (a1 < 1E-10) AND (delta > eql_TOLX) then msgbox("Bad guess(es)"); return({1,0}); end;
  if delta<eql_TOLX then
    if (2*n1)^(.5) < eql_TOLF then msgbox("Zero"); return({1,1}); end;
    if numvars<>numeq then
      costol:=abs(dot(vecgrad,vecD))/(2*n2);
      costol:=costol^(.5);
      if costol < eql_TOLLSQ then msgbox("Minimum"); return({1,2}); end;
    end;
  end;
  return({0,0});
end;

// ᴇ

Edit: for the sake of completeness, here is a sample EqLib Note file:

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"
    }
  }
}