How do I make a general solver for a system of equations (e.g. thermodynamic processes) with different inputs and outputs?

Question

Background

Let's say I have a system of equations (e.g. thermodynamics equations), where the "knowns" and "unknowns" are subject to change, and the system of equations can also change based on e.g. the type of thermodynamic process (isothermal, isobaric, isochoric, adiabatic).

Knowns and Unknowns Subject to Change

Take $PV=nRT$. Case 1: If I know $P$, $V$, $n$, and $R$, then $T\rightarrow\frac{PV}{nR}$. Case 2: I know $V$, $T$, $n$, $R$, then $P\rightarrow\frac{nRT}{V}$.

The "Hard-Coded" Solution

An easy solution is:

eqn = P V = n R T;
soln1 = Solve[eqn, T];
soln2 = Solve[eqn, P];

but this can become overwhelming with many input and output variables and especially if the systems of equations are also subject to change.

Question

How do I make make a general solver that takes a system of equations and whatever inputs are supplied (with units) and outputs the best attempt at a solution based on those inputs?

Some SE examples

I think this kind of approach is applicable to the following examples:

• Work done in Isobaric Process
• Answer to Comparison between isobaric, isothermal and adiabatic expansion
• Finding the Enthalpy of an Ideal Gas given internal energy
• Adiabatic proccess and Carnot cycle in a photon gas
• A simple thermodynamic question

Update

2020-09-19 I finally came across two related SE questions:

• Work with equation involving different independent variables [closed]
• Working with tables of equations in different units

Answer 1

General Solver

Define Function(s) to Retrieve System of Equations

Clear[P,V,n,R,T];
Rval=QuantityMagnitude@UnitConvert@Quantity[1, "MolarGasConstant"];
idealGasEqn := Module[{R=Rval,eqns}, eqns = {P*V == n*R*T}]

Known Variables

Case 1: P, V, and n are knowns (Solve for T)

Pval1 = Quantity[1.5, "Atmospheres"];
Vval1 = Quantity[3, "Liters"];
nval1 = Quantity[1, "Moles"];

Case 2: V, T, and n are knowns (Solve for P)

Vval2 = Quantity[3, "Liters"];
nval2 = Quantity[1, "Moles"];
Tval2 = Quantity[55,"Kelvins"];

Procedure

Setup

Equations, solve variables, and inputs

• Get system of equations based on input argument (e.g. type = "IdealGas") using a Switch statement.
• Define list of solve variables (Symbols that are left unset)
• Define list of input variables (mixture of set and unset)

Units

• Get output unit and SI unit Quantities both with magnitude 1
• Find positions in solve variable and input lists based on variable type (Symbol or Quantity) using Position
• Replace quantities with magnitude of SI-converted quantities using ReplacePart

Solve

Unitless Solution

• Solve for unknowns using SI magnitudes via UnitConvert and Solve

Unit Solution

• Attach SI magnitudes to unitless solution and convert to output units

Output

• Output a unitless or unit-containing solution

Module

idealGasSolver[P1_,V1_,n1_,T1_,type_:"IdealGas",unitlessQ_:False] := 
Module[
{eqns,vars},
(*get system of equations*)
eqns = Switch[type,"IdealGas",idealGasEqn];

vars = {P,V,n,T}; (*Symbols for solve, keep unassigned throughout*)
valsTmp = {P1,V1,n1,T1}; (*input values, some are Symbols, some are Quantities*)

(*units with magnitude 1*)
outUnits = Quantity[1,#]&/@{"Atmospheres","Liters","Moles","DegreesCelsius"};
SIunits = Quantity[1,#]&/@QuantityUnit@UnitConvert@outUnits;

(*find positions based on variable type*)
getIDs[head_] := Position[Head@#===head&/@valsTmp,True]//Flatten;
quantityIDs = getIDs[Quantity];
symbolIDs = getIDs[Symbol];

(*replace quantities with magnitude of SI - converted quantities*)
rules1 = MapThread[#1->#2&, 
 {quantityIDs,QuantityMagnitude@UnitConvert@valsTmp[[quantityIDs]]}];
vals = ReplacePart[valsTmp,rules1];

(*solve for unknowns using SI magnitudes, no units in output*)
rules2 = MapThread[#1->#2&,{vars[[quantityIDs]],vals[[quantityIDs]]}];
unitlessSoln = Solve[eqns/.rules2,vars[[symbolIDs]]][[1]];

(*convert solution to output units and include units*)
rules3 = MapThread[#1 -> #2 &, {vars[[symbolIDs]], 
 vals[[symbolIDs]]*SIunits[[symbolIDs]]}];
outVals = MapThread[UnitConvert[#1,#2]&, 
 {vars[[symbolIDs]]/.rules3/.unitlessSoln,outUnits[[symbolIDs]]}];
unitSoln = MapThread[#1->#2&,{vals[[symbolIDs]],outVals}];

(*output a solution based on unitlessQ argument*)
outsoln = If[unitlessQ,unitlessSoln,unitSoln]
]

Testing

Case 1

Clear[T];
idealGasSolver[Pval1, Vval1, nval1, T] (*output in units based on outUnits (deg C)*)
idealGasSolver[Pval1, Vval1, nval1, T, "IdealGas", True] (*output temperature SI unit (K) magnitude*)

{T -> Quantity[-218.31031631383098, "DegreesCelsius"]}
{T -> 54.83968368616898}

We get units with the first output, and an SI magnitude with the second.

Case 2

Clear[P];
idealGasSolver[P, Vval2, nval2, Tval2] (*output in units based on outUnits (atm)*)
idealGasSolver[P, Vval2, nval2, Tval2, "IdealGas", True] (*output pressure SI unit (Pa) magnitude*)

{P -> Quantity[2286477219992141/1519875000000000, "Atmospheres"]}
{P -> 2286477219992141/15000000000}

Exact arithmetic is preserved in this case.

Case 3 (additional case, underdetermined system of equations)

Clear[P, V, n]
idealGasSolver[P, V, n, Tval2] // N
idealGasSolver[P, V, n, Tval2, "IdealGas", True] // N (*output SI magnitude*)

{P -> UnitConvert[P*Quantity[1., "Kilograms"/("Meters"*"Seconds"^2)], Quantity[1., "Atmospheres"]], V -> UnitConvert[V*Quantity[1., "Meters"^3], Quantity[1., "Liters"]], n -> UnitConvert[P*V*Quantity[0.002186770091685928, "Moles"], Quantity[1., "Moles"]]}
{n -> 0.002186770091685928*P*V}

The second output (SI magnitude) is more parsable and less subject to issues if you were to apply this process successively (i.e. use the outputs as inputs to the next system of equations).