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

[Han]

Here's a very basic skeleton. This version does not handle systems that have free variables. Ideally, in the case of free variables, then the initial guesses should be used to establish a solution. For example, if (x,y,z) = (x,2x-1,3x+1), then the initial guess of (4,5,6) should produce the solution (4,7,13) using x=4. This will be implemented over time. For now, what you get is basically a proof of concept.

Code:

export eqldata;
export varlist;
export eqlist;
export eql_results;
export inputstr;

eql_init();
eql_cmenu();
eql_setup();
eql_dovars();
eql_dosolve();

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
    eqldata:=expr(tmp);
  end;
end;

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

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

  return(c);
end;

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

  for i from 1 to n do
    tmp:=tmp+"{eqlist(" + 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;

eql_dovars(c,f)
begin
  local eqns:=eqldata(c);
  local tmp:="";
  local n:=size(eqns(3));
  local i; local t;
  if f then
    varlist:=makelist(0,X,1,n);
  end;

  for i from 1 to n do
    tmp:="type(" + eqns(3,i) + ")";
    if CAS(tmp)==6 then tmp:=eqns(3,i) + ":=0"; CAS(tmp); end;
  end;

  tmp:="input({";

  for i from 1 to n do
    tmp:=tmp + "{" + eqns(3,i) + ",[0],{30,30," + string(i-1,1,12) + "}},";
    tmp:=tmp + "{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 + CAS(eqns(3,i)) + "," + 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 varlist(i)==0 then
        tmp:="purge(" + eqns(3,i) + ")";
        CAS(tmp);
      end;
    end;
  end;
  
  return(t);
end;

eql_dosolve(c)
begin
  local eqns:=eqldata(c);
  local n:=size(eqns(2));
  local i;
  local tmp:="fsolve([";
  local k:=ΣLIST(eqlist);

  for i from 1 to n do
    if eqlist(i) then
      k:=k-1;
      tmp:=tmp + eqns(2,i);
      if k then tmp:=tmp + ","; end;
    end;
  end;
  tmp:=tmp + "],[";

  n:=size(eqns(3));
  k:=n-ΣLIST(varlist);

  for i from 1 to n do
    if varlist(i)==0 then
      k:=k-1;
      tmp:=tmp + eqns(3,i);
      if k then tmp:=tmp + ","; end;
    end;
  end;
  tmp:="eql_results:=" + tmp + "])";
  inputstr:=tmp;
  CAS(tmp);
  local l:=size(eql_results);

  if type(eql_results)==6 then msgbox("No solution or infinite solutions; not implemented"); end;

// if multiple solutions; just get the first one; fix this later to choose "best" solution
  if type(eql_results)==4 then
    k:=1;
    for i from 1 to n do
      if varlist(i)==0 then
        tmp:=eqns(3,i) + ":=" + eql_results(1,k);
        CAS(tmp);
        k:=k+1;
      end;
    end;
  end;

end;

You must first create a Note named "EqLib" whose contents is:

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

Run eqlib() and choose "Linear Motion" (the only option). The next screen allows you to pick which equations to use. Leave all selected. The screen after that allows you to set up the known values. Those which should be treated as constants can be checked so that they are not considered a variable in the system. Press OK. If a solution exists, the values are updated. Press cancel to quit. I chose x=1, x0=5, t=3, and a=9.8. This results in v0=-16.033333... and v=13.366666...

Again, this is just a proof of concept. There are still a number of things I have to figure out about fsolve. Even when a system is solvable, if we give fsolve initial conditions, we sometimes get "Invalid dimension" errors. It appears that if fsolve is given more equations than the number of variables, it claims invalid dimension even though that may not be the case (e.g. some of the equations are equivalent).

An alternative to fsolve would be to implement our own system solver using something similar to that in SOLVESYS for the HP48/HP49. A brief glance at the source code suggests a Newton-like algorithm using the Jacobian.