systemSolver.frink

Download or view systemSolver.frink in plain text format


// Attempt to solve a system of equations.

use solvingTransformations.frink
showApproximations[false]
symbolicMode[true]

solveSystem[equations, solveForString] :=
{
   eqArray = new array
   size = length[equations]
   allUnknowns = new set
   for i=0 to size-1
   {
      syms = getSymbols[parseToExpression[equations@i]]
      equation = equations@i
//      println["Equation is $equation"]
//      println["Symbols are $syms"]
      eqArray.push[ [equation, syms] ]
//      println["eqArray is $eqArray"]
      allUnknowns = union[allUnknowns, syms]
//      println["allUnknowns is $allUnknowns"]
   }

//   println["All unknowns: $allUnknowns"]
   // Sort equations by number of unknowns (fewest first)
   sort[eqArray, {|a,b| length[a@1] <=> length[b@1]}]

   res = solveParts[eqArray, solveForString]
   res = eliminateDuplicates[res]
   res = eliminateOverconstrained[res]
//   res = eliminateSelfReferential[res]
   return res
}

// This is a recursive internal method to solve for other unknowns.
solveParts[eqArray, solveForString] :=
{
   results = new array
   size = length[eqArray]
   for i=0 to size-1
   {
      [eq, unknowns] = eqArray@i
      if unknowns.contains[solveForString]
      {
         oe = parseToExpression[eq]
         equation = parseToExpression["solve[$eq, $solveForString]"]
         solvedEq = transformExpression["solving", equation]
         results.push[solvedEq]
         println["Result is $solvedEq"]
         otherEqs = eqArray.shallowCopy[]
         otherEqs.remove[i]

         otherSize = length[otherEqs]
         for unknown = unknowns
         {
            pUnknown = parseToExpression[unknown]
            if (unknown != solveForString)
               for j=0 to otherSize-1
               {
                  res2 = solveParts[otherEqs, unknown]
//                  println["Unknown is $unknown"]
//                  println["Res2 is $res2"]
                  for respart = res2
                  {
//                     println["Replacing in $solvedEq, $pUnknown becomes " + child[respart,1]]
                    // TODO:  Replace solving with simplification rules.
                     r = transformExpression["solving",substituteExpression[solvedEq, pUnknown, child[respart,1]]]
//                     println["  Result: $r"]
                     results.push[r]
                  }
               }
         }
      }
   }
   return flatten[results]
}


eliminateDuplicates[eqArray] :=
{
   i=0
   j=1
   while (i<length[eqArray])
   {
      ie = eqArray@i
      while (j<length[eqArray])
      {
         if structureEquals[ie, eqArray@j]
         {
//            println["Removing duplicate " + eqArray@j]
            eqArray.remove[j]   // Don't advance index in this case
         } else
         j=j+1
      }
      i=i+1
      j=i+1
   }
   return eqArray
}

// This function eliminates overconstrained equations.  For example, a system
// containing the solutions a===1/2 c r  and  a===c d^-1 r^2  is
// overconstrained because a value can always be obtained with the first
// equation.  The second is not necessary.
eliminateOverconstrained[eqArray] :=
{
   size = length[eqArray]
   unknowns = new array
   for i = 0 to size-1
      unknowns@i = getSymbols[child[eqArray@i,1]]

   res = new array

   isProper = true
   ILOOP:
   for i=0 to size-1
   {
      overconstrained = false
      j = 0
      do
      {
         overconstrained = isProperSubset[unknowns@j, unknowns@i]
         if (overconstrained)
            println[eqArray@j + " is a proper subset of " + eqArray@i]
         j=j+1
      } while (j < size) && ! overconstrained

      if (! overconstrained)
         res.push[eqArray@i]  // If we got here, no j is a proper subset of i.
   }

   return res
}


symbolicMode[true]
//println[join["\n",solveSystem[["d === 2 r", "c === pi d", "a === pi r^2", "e===f", "g===h"], "r"]]]

// See https://web.archive.org/web/20210421140341/http://answers.yahoo.com/question/index?qid=20091120001614AAInec3
//
// Q. Given that P(A)=0.75, P(B|A)=0.8 and P(B|A')=0.6, calculate P(B) and P(A|B)
// Ans. 0.75, 0.8
// TODO:
//   eliminate pAprime and pAB, as those are what we want to solve for.
//println[join["\n",solveSystem[["pAB === pBA pA / pB", "pA === 1 - pAprime", "pAB === ( pBA pA ) / ( pBA pA + pBAprime pAprime ) "], "pB"]]]


Download or view systemSolver.frink in plain text format

This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 20145 days, 13 hours, 47 minutes ago.