What is the best way to parse differential equations and boundary conditions to a custom function?
Assuming everything is syntactically correct, the function that does what you want is Internal`ProcessEquations`SeparateEquations:
Internal`ProcessEquations`SeparateEquations[{y''[x] == y[x],
y[0] == 1, y[1] == 1}, {x}, {y}]
(*
{{}, {y[0] == 1, y[1] == 1}, {}, {y''[x] == y[x]}}
*)
Internal`ProcessEquations`SeparateEquations[
{y''[x] - z[x] == 0, y[0] == 1, y[1] == 1, z''[x] - y[x] == 0,
z[0] == 1, z[1] == 1}, {x}, {y, z}]
(*
{{}, {y[0] == 1, y[1] == 1, z[0] == 1, z[1] == 1},
{}, {-z[x] + y''[x] == 0, -y[x] + z''[x] == 0}}
*)
It's undocumented, and this appears to be its syntax and return value:
Internal`ProcessEquations`SeparateEquations[
{ equations },
{ indendent variables },
{ dependent variables }] (* N.B. No args: y not y[x] *)
(*
{ constraints on independent variables,
constraints on dependent variables, (* e.g BCs *)
algebraic equations,
differential equations }
*)
I've used this to write a parser that returns a data structure like @Nasser's. I don't mind sharing the code, but it is darn long, and I don't want to do too much refactoring to narrow its focus on your requirements.
Appendix: Parser Code Dump
The parser parseDiffEq[] is a somewhat pared-down version of the one alluded to above. It works with standard NDSolve input (omitting options):
myDE = parseDiffEq[{y''[x] == y[x], y[0] == 1, y[1] == 1},
y[x], {x, -1, 1}]
(*
<|"de" -> {y''[x] == y[x]},
"dependentVars" -> {y},
"independentVars" -> {x},
"completeDependentVars" -> {{y,y'}},
"bcs" -> {y[0] == 1, y[1] == 1},
"domain" -> {-1., 1.},
"return" -> y[x],
"firstorder" -> {y[1]'[x] == y[0][x], y[0]'[x] == y[1][x]},
"order" -> {{2}},
"type" -> "ODE"|>
*)
I cut some data structure items but left some that are not needed here but might be of interest.
The utility linearQ[], which will check if the DE is a linear system, seemed worth including given the OP's goal.
linearQ@myDE
(* True *)
Second example, a system:
my2DE = parseDiffEq[{y''[x] - z[x] == 0, y[0] == 1, y[1] == 1,
z''[x] - y[x] == 0, z[0] == 1, z[1] == 1}, {y[x], z[x]}, {x, -1, 1}]
(*
<|"de" -> {-z[x] + y''[x] == 0, -y[x] + z''[x] == 0},
"dependentVars" -> {y, z},
"independentVars" -> {x},
"completeDependentVars" -> {{y, y'}, {z, z'}},
"bcs" -> {y[0] == 1, y[1] == 1, z[0] == 1, z[1] == 1},
"domain" -> {-1., 1.},
"return" -> {y[x], z[x]},
"firstorder" -> {
-z[0][x] +y[1]'[x] == 0, -y[0][x] + z[1]'[x] == 0,
y[0]'[x] == y[1][x], z[0]'[x] == z[1][x]},
"order" -> {{2}, {2}},
"type" -> "ODE"|>
*)
linearQ@my2DE
(* True *)
Parser and utility code
There are internal, undocumented helper functions used that might be of interest:
Internal`ProcessEquations`SeparateEquations
Internal`ProcessEquations`FindDependentVariables
Internal`ProcessEquations`FirstOrderize
Internal`ProcessEquations`DifferentialOrder
Since they are undocumented, my ability to explain them is limited. The input to parseDiffEq[] is validated to some extend, but there are some checks I haven't gotten around to writing. The parser might occasionally fail on bad input without indicating why.
$parseKeys = { (* just a way for me to remember the data structure *)
"de", (* the diff. eqns. *)
"dependentVars", (* the "X" argument *)
"independentVars", (* the "Y" argument *)
"completeDependentVars", (* including lower-order derivatives *)
"bcs", (* boundary/initial conditions *)
"domain", (* interval of integration *)
"return", (* return expression *)
"firstorder",(* first-order equivalent system -- unnecessary *)
"order", (* differential orders of the DEs *)
"type" (* ODE, PDE,... -- unnecessary *)
};
ClearAll[parseDiffEq];
SetAttributes[parseDiffEq, HoldAll];
Options[parseDiffEq] = Thread[$parseKeys -> Automatic];
parseDiffEq::ndnl = NDSolve::ndnl;
parseDiffEq::dsvar = NDSolve::dsvar;
parseDiffEq::ndlim = NDSolve::ndlim;
(*
* Utilities
*)
ClearAll[
parseInterval, (* check indep var spec *)
validVariableQ, (* check whether an expression is a valid var *)
cullArgs, (* returns arguments of vars: y'[2]==0 -> {2} *)
varsToIndexedVars, (* convert Derivative[n][y] to y[n] *)
linearQ]; (* test whether a DE is linear *)
(* converts derivative y^(n) to y[n] *)
(* Used here for constructing the first order system
* and therefore unnecessary. Useful in other use cases
* for replacing derivatives by independent variables. *)
varsToIndexedVars[vars_][expr_] := varsToIndexedVars[expr, vars];
varsToIndexedVars[expr_, vars_] :=
With[{v = Alternatives @@ Flatten@{vars}},
expr /. {Derivative[n_][y : v] :> y[n], y : v :> y[0]}
];
(* taken from somewhere I've lost track of *)
validVariableQ[var_] := ! NumericQ[var] &&
FreeQ[var,
DirectedInfinity | Indeterminate] &&
(MemberQ[{Symbol, Subscript, K, C},
Head[var]] || ! AtomQ@Head[var] ||
Context[Evaluate@Head[var]] =!= "System`") &&
If[Head@Head[var] === Symbol,
! MemberQ[Attributes[Evaluate@Head[var]], NumericFunction],
validVariableQ[Head[var]]];
(* cullArgs - cull args of functions ff: {{args f1}, {args f2},..} *)
(* cullArgs[{y[0]==0,y[1]==0,z[0]==1},{y,z}] --> {{{0},{1}},{{0}}} *)
cullArgs[expr_, ff_] := DeleteDuplicates /@ Flatten[
Last@Reap[
Cases[
expr, (f : Alternatives @@ ff)[
args__] | _Derivative[f : Alternatives @@ ff][args__] :>
Sow[{args}, f], Infinity],
ff
],
1];
cullArgs[ff_][expr_] := cullArgs[expr, ff];
(* Checks if data structure de represents a linear equation or system *)
linearQ::usage = "linearQ[de] returns whether de is linear.";
linearQ[de_] := AllTrue[
Lookup[de, "de"],
Internal`LinearQ[
#,
Through[Flatten@{Lookup[de, "completeDependentVars"],
MapThread[
(Derivative @@ #2)@# &,
{Lookup[de, "dependentVars"], Lookup[de, "order"]}]} @@
Lookup[de, "independentVars"]]
] &
];
(* breaks down iterator {x,...} to {x, interval} and
* checks that x is a valid variable *)
parseInterval[xx : {x_, a___}] :=
If[! validVariableQ@x,
Message[parseDiffEq::dsvar, x];
Return[$Failed],
{x, {a}}
];
parseInterval[x_] := parseInterval@{x};
(*** end of utilities ***)
(*
* Main function: parses DE, vars, interval into an association
*
* Part I parses NDSolve style input into a sequence of option rules
* Part II construct the data struction Association[] from rules
*)
(* part I: parse equation and args into parts *)
parseDiffEq[eqns_List, yy_, xx_, deOpts : OptionsPattern[]] :=
Module[{
x, y, endpoints, interval,
conind, condep, alg, diff},
x = parseInterval@xx;
If[x =!= $Failed,
{x, interval} = x; (* split indep var and interval *)
y = yy /. v_[x] :> v; (*
strip arguments of dep var *)
{conind, condep, alg, diff} =
Internal`ProcessEquations`SeparateEquations[eqns, Flatten@{x},
Flatten@{y}];
(* TBD check validity {conind,condep,alg,diff} *)
endpoints = cullArgs[condep, Flatten@{y}];
interval = Flatten[{interval, endpoints}];
If[Length@interval == 0,
Message[parseDiffEq::ndlim, xx];
x = $Failed,
If[! VectorQ[interval, NumericQ],
Message[parseDiffEq::ndnl,
First@Cases[interval, x0_?(! NumericQ[#] &)], interval];
x = $Failed,
interval = MinMax@N@interval (* N[] optional;
use WorkingPrecision? *)
]
]
];
parseDiffEq[
"de" -> diff,
"bcs" -> (condep /. Automatic -> {}),
"independentVars" -> Flatten@{x},
"dependentVars" -> Flatten@{y},
"return" -> yy,
"domain" -> interval,
deOpts] /; FreeQ[x, $Failed]
];
(* part II: check and process parts given as option rules *)
parseDiffEq[opts : OptionsPattern[]] :=
Module[{asc, alldvars, firstordersys, foRules},
(* TBD: validate option values ??? *)
(** set up association from options **)
asc = <|Thread[$parseKeys -> OptionValue@$parseKeys]|>;
(** parses indep var from eqns; NDSolve does not do this -- unnecessary **)
If[asc@"independentVars" === Automatic,
asc@"independentVars" =
DeleteDuplicates@
Cases[Flatten@{asc@"de"}, _[x__Symbol] |
Derivative[__][_][x__Symbol] :> x, Infinity]
];
(** check type of DE -- unnecessary **)
asc@"type" = Switch[Length@asc@"independentVars"
, 0, "Algebraic" (* unsupported *)
, 1, "ODE"
, n_Integer /; n > 1, "PDE" (* unsupported *)
, _, $Failed];
(** parse dependend variables from equations -- unnecesary **)
If[asc@"dependentVars" === Automatic
, asc@"dependentVars" =
Internal`ProcessEquations`FindDependentVariables[
Flatten@{asc@"de"}, asc@"independentVars"]
];
(** construct first-order equivalent system -- unnecessary **)
firstordersys =
Internal`ProcessEquations`FirstOrderize[#1, #2, 1, #3] & @@
Lookup[asc, {"de", "independentVars", "dependentVars"}];
alldvars = firstordersys[[3]] /. firstordersys[[4]];
If[VectorQ[alldvars], alldvars = List /@ alldvars];
asc@"completeDependentVars" = alldvars;
foRules =
MapAt[ (* replaces NDSolve`y$nnn$1 by y[1] etc *)
varsToIndexedVars[Lookup[asc, "dependentVars"]],
Flatten@{firstordersys[[4]], # -> # & /@
Lookup[asc, "dependentVars"]},
{All, 2}];
asc@"firstorder" =
Join[firstordersys[[1]], firstordersys[[2]]] /. foRules;
(** store differential order -- unnecessary **)
asc@"order" =
Internal`ProcessEquations`DifferentialOrder @@
Lookup[asc, {"de", "independentVars", "dependentVars"}];
asc
];
I'll just give an idea to make it easier to do this. Which is not to use the same API as NDSolve, as that requires much much more work to parse it.
Instead, have the caller pass the input in Association.
Yes, this might be a little bit more work for the user, but not much. On the other hand, this greatly simplifies the parsing and checking inside your ndsolve, because now all the entries can be accessed directly by field names from the association instead of using pattern search.
This is actually how number of other software do it. The user fills in a "record" or a "struct" in C talk, and passes this struct to the function to process.
The function then just reads the values directly from the record by name.
There is a quick prototype. This will work for many number of odes.
You build one association for each ode
ClearAll[y, x, z, ode1, ode2];
ode1 = <|"depVar" -> y,
"indepVar" -> x,
"ode" -> y''[x] - z[x] == 0,
"ic" -> {y[0] == 1, y[1] == 1}|>;
ode2 = <|"depVar" -> z,
"indepVar" -> x,
"ode" -> z''[x] - y[x] == 0,
"ic" -> {z[0] == 1, z[1] == 1}|>;
domain = {{x, -1, 1}};
setOfODES = {ode1, ode2};
Now you call your ndsolve
ndsolve[setOfODES, domain]
And this is ndsolve
ndsolve[odes_List, domain_List] := Module[{n = Length@odes, m, currentODE},
Print["You entered ", n, " odes"];
Do[
currentODE = odes[[m]];
Print["\nODE ", m, " is ", currentODE["ode"],
"\nthe dependent variable is ", currentODE["depVar"],
"\nthe independent variable is ", currentODE["indepVar"]
]
, {m, 1, n}
];
(*example how to read fields from association*)
If[n > 1,
If[ Length@Union["indepVar" /. odes] > 1,
Return["Error, independent variable must be the same", Module]
]
];
(*many many more additional checks and balances*)
(*check domain is valid*)
(*check initial conditions are valid and using same symbols,etc...*)
Print["Parsed OK"]
(*now you can go on and actually numerically solve them. But the hard work*)
(*has been done above, which is parsing, the rest is easy :) *)
]
And it gives this output
You entered 2 odes
ODE 1 is -z[x]+y''[x]==0
the dependent variable is y
the independent variable is x
ODE 2 is -y[x]+z''[x]==0
the dependent variable is z
the independent variable is x
Parsed OK
The above is just the start. But the main point, it is much easier now to handle, since you do not have to do too much parsing, compared to the way NDSolve takes its input as lists, where you'd have to parse the content of each list, pick which part is which, and so on. This is at the cost, the caller has to set up an association for each ODE. But I think it is not a big deal to do.
Here's a simpler way (simpler than my first answer) that I came up with today exploring a problem with DSolve. It calls DSolveValue and intercepts the DSolve parser and returns an association with the equations broken down by type, before the system is solved:
parseODE@NDSolve[{y''[x] == y[x], y[0] == 1, y[1] == 1}, y[x], {x, -1, 1}]
(*
<|"OtherEquations" -> {}, (* nonempty => error (probably) *)
"BoundaryConditions" -> {y[0] == 1, y[1] == 1},
"Algebraic" -> {}, (* algebraic equations in terms of y and x *)
"Differential" -> {y''[x] == y[x]},
"Dependent" -> {y},
"Independent" -> {x},
"Region" -> {x, -1, 1}|> (* see the PDE example below *)
*)
Code for function:
ClearAll[parseODE];
SetAttributes[parseODE, HoldFirst];
$dsolvers =
DSolve | DSolveValue | NDSolve | NDSolveValue | ParametricNDSolve |
ParametricNDSolveValue;
parseODE[
_?(MatchQ[$dsolvers])[
eqns_, v_, t : Longest[{_, _?NumericQ, _?NumericQ} ..] | _, ___]
] := parseODE[eqns, v, t];
parseODE[eqns_, v_, t__] :=
Block[{DSolve`DSolveParser =
Function[{eqs, dependent, independent, stuff},
Return[
With[{independents =
Flatten@{independent /.
{{x_, _?NumericQ, _?
NumericQ} :> x, vv_ \[Element] _ :> vv}
}},
Join[
AssociationThread[
{"Other", "Initial", "Algebraic", "Differential"} ->
Internal`ProcessEquations`SeparateEquations[
Flatten@eqs, independents, dependent]],
<|"Dependent" -> dependent,
"Independent" -> independents,
"Region" -> independent|>
]],
Block]
]},
DSolveValue[eqns, v, t]
]
More examples. Note that the domain {x, 0, 1}, {t, 0, 1} for the PDE in the first example is rewritten by DSolveValue into an ImplicitRegion. The others show variation in input type (x instead of {x, 0, 1}, a system instead of a single ODE).
weqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}];
ic = {u[x, 0] == E^(-x^2), Derivative[0, 1][u][x, 0] == 1};
parseODE@DSolveValue[{weqn, ic}, u[x, t], {x, 0, 1}, {t, 0, 1}]
(*
<|"OtherEquations" -> {},
"BoundaryConditions" -> {{u[x, 0] == E^-x^2, Derivative[0, 1][u][x, 0] == 1}},
"Algebraic" -> {},
"Differential" -> {Derivative[0, 2][u][x, t] == Derivative[2, 0][u][x, t]},
"Dependent" -> {u},
"Independent" -> {x, t},
"Region" -> {{x, t} \[Element]
ImplicitRegion[0 <= x <= 1 && 0 <= t <= 1, {x, t}]}|>
*)
parseODE@DSolve[{y''[x] == y[x], y[0] == 1, y[1] == 1}, y[x], x]
(*
<|"OtherEquations" -> {},
"BoundaryConditions" -> {y[0] == 1, y[1] == 1}, "Algebraic" -> {},
"Differential" -> {y''[x] == y[x]},
"Dependent" -> {y}, "Independent" -> {x}, "Region" -> {x}|>
*)
parseODE@NDSolveValue[{a'[t] == 1, y'[t] == 1, a[0] == 0,
y[0] == 0}, {a[t], y[t]}, {t, 0, 1}]
(*
<|"OtherEquations" -> {},
"BoundaryConditions" -> {a[0] == 0, y[0] == 0}, "Algebraic" -> {},
"Differential" -> {Derivative[1][a][t] == 1,
Derivative[1][y][t] == 1}, "Dependent" -> {a, y},
"Independent" -> {t}, "Region" -> {t, 0, 1}|>
*)
If the differential order(s) of the variables would be useful, one could add a line to the association:
"Order" -> Internal`ProcessEquations`DifferentialOrder[
Flatten@eqs, independents, dependent]