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

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]

Answer 4

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.

Using

WolframAlpha["y''[x] == y[x],y[0] == 1, y[1] == 1"]

is simply the best. The advantages are superior to everthing else on the market.

It name the autonomous equation. It classifies the ODE. It solve the differential equation. You get option to calculate an approximate form and get a step by step solution. Depending on the importance of the ODE there is more information. Even a plot and some variation on solutions parameters can be presented. This is traditional mathematica notation and therefore more common. The solution of the parsing address much more audiance than Mathematica code. Using the pods as I introduced some of them does much more than all the other answerers to the question do. The knowledge is much more rich and versatile.

It depends on the interpretation of parsing. My answer follows the ideation and concepts of Wolfram Alpha complete. Not other solution does this. And the other answer are not parsing in the sense of using knowledge database. They use information otherwise hidden by the built-ins because they are tautologies derive solely from the input by the user.

Make use for example of the definition from Parsing. With my answer the collection of answers to the question start to match what parsing is. This is the step from simple computational language to knowledge databased natural language parsing. The step from 4GL to 4GL knowledge based natural parsing.

The point can be that the other use backtracking of the input, while my solution is real parsing of tokens in natural language mapped on a rich knowledge database.

My solution offers access to the step-by-step-solution of Wolfram Alpha while all the others only envolve the Notebook-In-Out-schema.

WolframAlpha["y''[x] - z[x] == 0, z''[x] - y[x] == 0"]

WolframAlpha give the ODE classification system of ordinary differential equation. Indeed this is second order and linear too. This classification is interpreted the rest is automated knowledge retrieval from the knowledge database about differential equations.

Since none of the is probably new a solution is always a reproduction of some to human mankind known problem of ordinary differential equations. The reproduction of solutions may either be done by algorithms rearranging in the input or graph based search algorithms.

The internal built-ins

Internal`ProcessEquations`SeparateEquations
Internal`ProcessEquations`FindDependentVariables
Internal`ProcessEquations`FirstOrderize
Internal`ProcessEquations`DifferentialOrder

are examples of such categories of graph keys to organize knowledge in need a priori for a efficient and correct solution. They are introduced earlier but reorganized in the knowledge methodology with the introduction into Mathematica. Both are not only compatible they rely on the very same methodologies internally.

It should be by the degree of maturity to Mathematica that Association are not very fast but most optimal for knowledge representation. Both represenation can be transformed into each without loss of information.