# HP Forums

Full Version: Creating an equation library (Updated 03-FEB-2017)
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Pages: 1 2 3
Hello,
Not an answer to why complex elements in the first matrix: just observed that the second vector element is nearly zero, so the complex elements have no impact on L1(1)*diag(L1(2))*TRN(L1(3)).
Seems like a floating point rounding issue in the imaginary part calculation of Σ which translates to some really bad error in U, so that it isn't unitary anymore (L1(1)*trn(L1(1))=trn(L1(1))*L1(1)=identity doesn't hold).
(01-15-2015 12:50 AM)Snorre Wrote: [ -> ]Hello,
Not an answer to why complex elements in the first matrix: just observed that the second vector element is nearly zero, so the complex elements have no impact on L1(1)*diag(L1(2))*TRN(L1(3)).
Seems like a floating point rounding issue in the imaginary part calculation of Σ which translates to some really bad error in U, so that it isn't unitary anymore (L1(1)*trn(L1(1))=trn(L1(1))*L1(1)=identity doesn't hold).

Maybe the algorithm uses 1/s_i * A * v_i where s_i is the i-th eigenvalue and v_i is the corresponding eigenvector (perhaps this is why the corresponding column of U is all complex -- but I can't seem to make the results match this...)

Another way to get non-diagonal matrices is to compute the orthonormal eigenvectors of A^T * A and of A * A^T, for which the Prime does not seem to introduce complex values. The problem is that for non-square matrices additional steps may be needed to ensure that the eigenvectors are orthogonal (usually Gram-Shmidt takes care of this).

Still, that's a lot of additional work that should have been automatically handled by SVD... :-(
Hello Han,
by making your matrix containing only reals (just change the very first element from 3 to 3.) you won't get complex results.
This new result is real but the first matrix is not orthogonal again (maybe a minor bug or -- worse -- just a numerically instable svd).
(01-15-2015 02:55 AM)Snorre Wrote: [ -> ]Hello Han,
by making your matrix containing only reals (just change the very first element from 3 to 3.) you won't get complex results.
This new result is real but the first matrix is not orthogonal again (maybe a minor bug or -- worse -- just a numerically instable svd).

I'm unable to make this happen; changing 3 to 3. still shows up as 3 in Home view; in CAS view it still gives me the same complex results (as well as Home). I use standard notation for the number format, if that makes any difference.
In GIAC/XCAS: evalf(SVD([[3,4.5,4.5],[4.5,12.25,1.25],[4.5,1.25,12.25]]))

no problems with complex elements in the U matrix, but the V matrix has all zeros in column 3 (in Wolfram alpha, for this example U equals V). The Sigma matrix, containing the singular values, is the same.
(01-15-2015 04:25 AM)Helge Gabert Wrote: [ -> ]In GIAC/XCAS: evalf(SVD([[3,4.5,4.5],[4.5,12.25,1.25],[4.5,1.25,12.25]]))

no problems with complex elements in the U matrix, but the V matrix has all zeros in column 3 (in Wolfram alpha, for this example U equals V). The Sigma matrix, containing the singular values, is the same.

Hm... that's clearly a bug in giac/xcas; for symmetric matrices the SVD should be essentially the same as an eigen-decomposition (up to a few nuances with signs of eigenvalues and corresponding eigenvectors).
Hello Han,
(01-15-2015 04:20 AM)Han Wrote: [ -> ]I'm unable to make this happen; changing 3 to 3. still shows up as 3 in Home view; in CAS view it still gives me the same complex results (as well as Home). I use standard notation for the number format, if that makes any difference.
Forgot to mention: that's possible by entering the matrix in CAS or do a approx([[...]]).
But you're right. It doesn't make any difference.

I did a mistake (luckily didn't clear the history): I've introduced some invisible roundoff deviation by
[attachment=1448]
Looks like the SVD command uses EIGENVV to obtain the matrix V and to determine the diagonal matrix of eigenvalues. Letting M0 be the matrix I mentioned earlier, then EIGENVV(M0*trn(M0)) gives

Code:
M0:=[[3,4.5,4.5],[4.5,12.25,1.25],[4.5,1.25,12.25]]; EIGENVV(M0*trn(M0));

The result:
Code:
{ [  [−0.426401432711,0.904534033733,4.07681407754ᴇ−17],  [−0.639602149067,−0.301511344578,0.707106781187],  [−0.639602149067,−0.301511344578,−0.707106781187] ], [  [272.25,0,0],  [0,−7.5863407865ᴇ−32,0],  [0,0,121] ] }

So the middle eigenvalue for M0*trn(M0) is a negative value. Take the square root and you get the eigenvalue found from doing SVD(M0). Note that the real part was dropped (likely due to roundoff).

(−7.5863407865E−32)^.5 = 1.48466797577e−30+2.75433127755e−16*i

I wonder if within the SVD algorithm, forcing the diagonal matrix from the eigen-decomposition to have all positive entries (they necessarily have to be, since they are supposed to be squares), would provide a quick fix for this issue.

There is, of course, always the option of writing our own SVD...
I hope the SVD algorithm will be checked by Bernard Parisse and HP, especially since it appears that XCAS is affected.

(Of course, SVD, being numerical in nature, is designed to be used in Home, so maybe not).
(01-16-2015 04:30 PM)Helge Gabert Wrote: [ -> ]I hope the SVD algorithm will be checked by Bernard Parisse and HP, especially since it appears that XCAS is affected.

(Of course, SVD, being numerical in nature, is designed to be used in Home, so maybe not).

I have a sneaky suspicion it will be fixed :-).
To create this eq library, is required a function that solves nonsquare systems (more unknowns than equations) to be devised an algorithm that recursively solving I think I saw one HPCalc
(01-18-2015 08:07 PM)compsystems Wrote: [ -> ]To create this eq library, is required a function that solves nonsquare systems (more unknowns than equations) to be devised an algorithm that recursively solving I think I saw one HPCalc

The plan is to use Newton's method (which is iterative). Each set of equations is transformed into functions $$\vec{F} = \langle f_1, f_2, \dotsm, f_m \rangle^T$$. Applying Newton's method to solve $$\vec{F}(\vec{x}) = 0$$ would require the Jacobian. Snorre posted an example in this thread. In order to solve for $$\vec{x}$$, we use SVD. This will take care of all cases, but one would have to properly interpret the solutions. SOLVESYS uses a similar scheme, but does not seem to use SVD or even LSQ for overdetermined systems.
Another option is to create a version SVD, I extracted the source code of (SVD MathTools TI 89) With tiedit http://home.arcor.de/p-engels/tiedit/
/!\
TM = TRANSPOSE
Code:
SVD(mm) Func ©SVD(mat) computes the singular value decomposition of numeric matrix mat ©Bhuvanesh Bhatt Local mm2,eps,evl,evc,isxact,nzevl,tmp,tmp2,tmp3,ss,uuu,vvv dim(mm)→tmp2 inString(string(mm),".")=0→isxact mm TM*mm→mm2 If not isxact Then   real(eigVl(mm2))→evl   ©"If dim(evl)>max(tmp2,tmp2) Then:mm*(mm)TM→mm2:real(eigVl(mm2))→evl:EndIf"   when(getType(tol)="NUM",tol,5.E­13,5.E­13)→eps   ©"seq(when(real(evl[ii])<eps,0,real(evl[ii]),real(evl[ii]))+when(imag(evl[ii])<eps,0,imag(evl[ii]),imag(evl[ii]))*tm,ii,1,dim(evl))→evl"   conj(eigVc(mm2))→evc Else   mathtool\eigenval(mm2)→evl   0→eps   conj(mathtool\gschmidt(mathtool\eigenvec(mm2)))TM→evc EndIf Σ(when(abs(evl[ii])>eps,1,0,1),ii,1,dim(evl))→nzevl mathtool\reverse(mathtool\sort(seq(augment({evl[ii]},mat▸list(evc[ii])),ii,1,dim(evl))))→tmp subMat(tmp,1,2)→vvv seq(tmp[ii,1],ii,1,rowDim(tmp))→evl √(left(evl,min(tmp2,tmp2)))→tmp3 diag(tmp3)→ss seq(mat▸list(tmp3[ii]^(­1)*mm*conj(vvv[ii])),ii,1,nzevl)→uuu If tmp2>nzevl   augment(uuu;subMat(identity(colDim(uuu)),1,1,tmp2-nzevl,colDim(uuu)))→uuu conj(mathtool\gschmidt(uuu TM))→uuu If tmp2>tmp2 Then   seq(seq(0,jj,1,tmp2),ii,1,tmp2-tmp2)→tmp3   If getType(tmp3)="LIST"     {tmp3}→tmp3   augment(ss;tmp3)→ss EndIf If tmp2<tmp2 Then   seq(seq(0,jj,1,tmp2-tmp2),ii,1,tmp2)→tmp3   If getType(tmp3)="LIST"     {tmp3}→tmp3   augment(ss,tmp3)→ss EndIf ©if rowdim(vvv)≠coldim(vvv) or rowdim(vvv)≠tmp2 or rowdim(uuu)≠coldim(uuu) or rowdim(uuu)≠tmp2 or dim(ss)≠tmp2 return "Error: Final matrix structure inconsistent. Please report to the author." {uu=uuu,ΣΣ=ss,vv=conj(vvv)} EndFunc

MathTools Guide

http://www.ibiblio.org/technicalc/packag...hTools.pdf

/////////

PD: A help from the community would be create a form to add equations, as doing FORMULAPRO for TInspire calculators

http://tiplanet.org/forum/viewtopic.php?t=9576 http://education.bwns.be/FormulaPro/
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:=; oldvecF:=; vecX:=; oldvecX:=; vecD:=; 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) + "),,{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:=;   vecX(numvars):=0; oldvecX:=vecX; vecD:=vecX;   vecF:=;   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"     }   } }
The next version will likely make use of the Solve app and that the end result would be a mix of an equation library (a la HP48) and SolveSys (for both HP48 and HP49). The choice of basing the next version on the Solve app is that the Symb view would allow one to view all the equations in textbook form. That said, there are some caveats with this idea:

1. Variables would have to be created first. This can be done easily as the user would simply be asked to enter a list of variable names, and their creation would be done by the app. While it is possible to simply create equations without specifying the variable names (and have the app create them where necessary), we would be limited to algebraic input and display of said equations. (In the Symb view, one cannot enter an expression/equation involving non-existent variables). For equations that involve complex expressions, I find that being able to view them using textbook display to be more advantageous than not having to specify the variables of a system.

2. We are limited to only 10 expressions/equations using the Symb view. This limitation is only in the viewing of the equations. The system would be stored internally by the app, and there would be mechanisms to "page" through all the equations in a system having more than 10 equations. I have not run into too many instances of systems involving more than 10 equations except in the case of regressional analysis (curve fitting data to a formula). In the latter case, all the equations would look essentially the same, though.

Once a user sets up the variables and equations, there will be options to the "save" the system and add it into one's "library" of equations. Then, users can either create a new system, or select from their library.

If anyone has other ideas they'd like to propose, please let me know.
Looks great!

One comment: Will there be an option to solve for a variable in C vs. restricting it to R (e.g., with a checkbox [Cmplx On/Off], kind of like SolveSys)?
(11-13-2015 08:05 PM)Helge Gabert Wrote: [ -> ]Looks great!

One comment: Will there be an option to solve for a variable in C vs. restricting it to R (e.g., with a checkbox [Cmplx On/Off], kind of like SolveSys)?

Yes, even the current solver handles complex solutions (or at least it should). If it doesn't then there's probably some place where dot products aren't using conjugates properly. CAS dot product is different from Home dot products.
Thanks - - it seems to work now - - still have to test it a little bit more.
(11-13-2015 03:20 PM)Han Wrote: [ -> ]The next version will likely make use of the Solve app and that the end result would be a mix of an equation library (a la HP48) and SolveSys (for both HP48 and HP49). The choice of basing the next version on the Solve app is that the Symb view would allow one to view all the equations in textbook form. That said, there are some caveats with this idea:

1. Variables would have to be created first. This can be done easily as the user would simply be asked to enter a list of variable names, and their creation would be done by the app. While it is possible to simply create equations without specifying the variable names (and have the app create them where necessary), we would be limited to algebraic input and display of said equations. (In the Symb view, one cannot enter an expression/equation involving non-existent variables). For equations that involve complex expressions, I find that being able to view them using textbook display to be more advantageous than not having to specify the variables of a system.

2. We are limited to only 10 expressions/equations using the Symb view. This limitation is only in the viewing of the equations. The system would be stored internally by the app, and there would be mechanisms to "page" through all the equations in a system having more than 10 equations. I have not run into too many instances of systems involving more than 10 equations except in the case of regressional analysis (curve fitting data to a formula). In the latter case, all the equations would look essentially the same, though.

Once a user sets up the variables and equations, there will be options to the "save" the system and add it into one's "library" of equations. Then, users can either create a new system, or select from their library.

If anyone has other ideas they'd like to propose, please let me know.

As I was reading your suggestions it seemd you were reading on my mind. Your targets are perfect, a hybrid of fast inserting or pre-defined collection of equations. Wonderful. This program would be considered one of my top- rated Hp program ever (with Vigag, FEM and Secciones). I´ll stay tuned for your release. Thank you Han for your great work.
So, in order to obtain the same DOT() result in CAS as in Home,

in CAS only,

one has to call DOT(CONJ([vector_1]),[vector_2]);

with complex elements in vector_1, in order to negate the complex conjugates which otherwise arise. A bug?
Pages: 1 2 3
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :