# 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:

# Update

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

# 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
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*)
quantityIDs = getIDs[Quantity];
symbolIDs = getIDs[Symbol];

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

(*solve for unknowns using SI magnitudes, no units in output*)
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]]}];
{vars[[symbolIDs]]/.rules3/.unitlessSoln,outUnits[[symbolIDs]]}];

(*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).