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

, 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:

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;

// ᴇ