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;
// ᴇ