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
];
how can I reproduce this (smart) behaviour in my home made function
by parsing the API. These things are not easy. There are two stages to parsing. First you have to make sure the input is valid. Then you have to pick the pieces you want from the input so they can be used internally. First stage is basically verification stage. It checks the input is all valid (right types, no conflicts, etc...). Once this is done, then parser collects all the pieces andsend them in right order to be processed by NDSolve (may be in different structure) So I am not sure what is it you asking to obtain here. $\endgroup$ – Nasser May 14 '20 at 9:25NDSolve::ndnco: The number of constraints (3) (initial conditions) is not equal to the total differential order of the system plus the number of discrete variables (4).
$\endgroup$ – Nasser May 14 '20 at 9:38