maxima CAS library for simplifying list of equations

applyN.mac

 /* put this in the same directory as the wx maxima file (.wxmx)
 load("applyN.mac");
 then use applyN as below
 */
 /*  applyN eliminates variables from a set of equations,
  *  while maintaining expressions for the eliminated variables.
  *
  * For example:
  *
  * r : applyN([v1 + v2 = v4^2], [v1])
  *
  * r[1] will be [v1 = v4^2 - v2]
  *
  * r[2] will be [], but if there were more equations they would be left there for
  * your inspection so that you can decide which variables can still be eliminated.
  *
  * Often an engineering model will have a large number of variables,
  * but the core of the model is a small system of nonlinear equations
  * with just a few degrees of freedom.
  * It is useful to keep expressions for the intermediate values, while also
  * reducing it to a small number of variables / constraint equations
  * before the next step which is either a numerical method like maxima's fmin_cobyla,
  * or plotting contours etc. with R ggplot
  */
applyN(eqs, vs) := lreduce(apply1, vs, [[], eqs])$

/* count the number of tokens */
ntokens(expr) := if atom(expr) then 1 else 1 + apply("+",map(ntokens,args(expr)));

goodsol(expr) := not(member(lhs(expr), listofvars(rhs(expr))));


deleteix(ix, xs) := append(firstn(xs, ix-1), rest(xs, ix))$

/* try to eliminate Bi (ie. apply1(eqs, x) will try to eliminate x) from eqs=[ks,es]
  *
  * this means solving each e in es for Bi,
  * choosing the simplest solution "s". Then
  * return the new [s: ks, es \\ e],
  * where the the solved equation is moved into the ks,
  * and all occurrences of Bi are removed in RHSs (by
  * substituting the solution)
  * */
apply1(eqs, Bi) := block([bi12, bi, bis, nonbi, subbi, iden, rhseqs],
        bis() := sublist(
                        create_list( block([s : linsolve([eqs[2][i]], [Bi])],
                                          if s = [] then false else [i, s[1]]),
                                        i, 1, length(eqs[2])),
                      lambda([x], x # false)),
        bi12 : first(sort(bis(), lambda([x, y], goodsol(x[2]) and (ntokens(x[2]) < ntokens(y[2]))))),

        bi : bi12[2],

        subbi : lambda([v], subst(bi, v)),
        /* nonbi is wrong because lhs(b) # Bi
          is only right if assoc(Bi, eqs[2]) is okay,
          I need to propagate the index into eqs[2]
          that led to the solution selected...
          I think I fixed that but it still isn't right
          Now I suspect that the problem is that I am
          getting a bi with whose rhs contains the lhs,
          and so the number of variables isn't reduced while
          the number of equations is.

          sublist is wrong I want to delete the first one

          for now the workaround is to adjust eqs[2] is to replace
          [f = x + y, f = y + z] with something like
          [f = x + y, f - z = y]
        */
        nonbi : if member(Bi, listofvars(rhs(bi))) then lambda([v], v) else
                if listp(bi12) then lambda([v],deleteix(bi12[1], v)) else
                lambda([v],sublist(v, lambda([x], lhs(x) # Bi))),
        iden : sublist_indices(eqs[2], lambda([x], lhs(x) = Bi)),
        rhseqs : if length(iden) > 1 then create_list(iden[i] - iden[i+1], i, 1, length(iden)-1) else [],
            /* if I have   [ Bi = x, Bi = y, Bi = z ],
               bi will Bi = x, though perhaps linsolve should always be tried
               and nonbi(eqs[2]) will be empty,
               which loses two equations  x = y, y = z that still can be considered
                */
        [cons(bi, subbi(eqs[1])),
         subbi(append(rhseqs,
          nonbi(
                 eqs[2])))]
        )$

/* additional definitions for fortran export */
optimprefix : e$

/* maxima equivalent of Data.List.\\ */
listdiff(a,b) := lreduce(lambda([x_, y], delete(y, x_)), b, a)$

/* maxima equivalent of Data.List.intersperse */
intersperse(_, L) := if emptyp(L) then []
  elseif emptyp(rest(L)) then [first(L)]
  else append([first(L), _], intersperse(_, rest(L)));

sprintlns([exprs]) :=  apply(sprint, intersperse(newline, exprs))$

fortranHeader(name, invars, outvars) := block([name1 : concat(name, ".f90")],
  with_stdout(name1, sprintlns(concat("subroutine ", name, "(vars, result)"),
  "  implicit real*8(a-z)",
  concat("  real*8, intent(in) vars(" , length(vars) , ")"),
  concat("  real*8, intent(out) result(",length(outvars),  ")"),
  concat("  INCLUDE ", name, ".inc"),
  concat("end subroutine ", name))
  ))$