What is the best way to parse differential equations and boundary conditions to a custom function?

Question

I have written a method to turn systems of linear differential equations into matrix equations (discretisation). This handles boundary conditions using the row replacement method.

At the moment, I am trying to find the best way to build this method into a custom NDSolve type of callable function with both equations and boundary conditions provided as arguments.

For concreteness, imagine I want to solve $y''(x)=y(x)$ on the interval $[-1,1]$ with conditions $y(0)=1$ and $y(1)=1$ (notice that the first condition is not at the boundary). If I were to use NDSolve, I would simply input:

NDSolve[{y''[x] == y[x], y[0] == 1, y[1] == 1}, y[x], {x, -1, 1}]

And NDSolve would immediately interpret the first element of the list as the differential equation, and the two others as constraints ("boundary" conditions).

For a system of two equations:

NDSolve[{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}]

would work equally well.

My question is thus : how can I reproduce this (smart) behaviour in my home made function ? How does NDSolve handle the parsing of arguments ?

Answer 1

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
   ];

Answer 2

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.

Answer 3

The best available off-the-shelf solution for this task shall be:

WolframAlpha["y'[x]\[Equal]x^3", IncludePods -> "ODEClassification", 
 AppearanceElements -> {"Pods"}]

__

This can be driven further:

{WolframAlpha["y'[x]\[Equal]x^3", {{"ODENames", 1}, "Content"}], 
 WolframAlpha["y'[x]\[Equal]x^3", {{"ODENames", 2}, "Content"}]}

WolframAlpha["y'[x]\[Equal]x^3", IncludePods -> "ODENames", 
 AppearanceElements -> {"Pods"}]

---

It is somehow a harsh reduction to life with names because this is a type of ODE and applicable to a lot of equation. Names would have been nice for the typing

G[y',x]==0

with

G=y'-x^3

Or both functions in y' and x^3 are polynomials of positive integer order.

The step-by-step solution in Mathematica is somewhat too short.

x^3dx

integrates in general to (this might need the Integrate built-in)

x^4/4+constant

dy

integrates in general to (this might need the Integrate built-in)

y

So in the equation

y[x_]:=x^4/4+constant

Since exact and separable are the more important typings, attribute of the given ODE they have to be favoured over

first-order linear ordinary differential equation

Exact implies there is a closed solution that is critical case is based on inverse functions. That does not follow from this type classification.

There is a need to discuss the symmetry of the two function. In this case both odd, so the resulting solution functions will be even.

If a comparison is initiated with DETools from Maple there are many question open:

info = WolframAlpha[
   "(x^2+x Exp[x])D[y[x],x,x,x,x]+(4 Exp[x]+6 x+ 3 Exp[x] \
x)D[y[x],x,x,x]+(6 + 3 Exp[x] x + 9 Exp[x])D[y[x],x,x]+(Exp[x] x+ 6 \
Exp[x])D[y[x],x]+Exp[x]y[x]\[Equal]0", "PodInformation"];
ids = Rest[DeleteDuplicates[#[[1, 1, 1]] & /@ info]];
titles = Map[{{#, 0}, "Title"} &, ids] /. info;
contents = 
  Column[Cases[info, _[{{#, _}, "Content"}, val_] :> val]] & /@ ids;
MenuView[Thread[titles -> contents], ImageSize -> Automatic]

This is following odeadvisor a fully exact ordinary differential equation, linear and higher-order. Mathematica avoids the pod with the ODE names. There are more differences.

stj