3
$\begingroup$

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 ?

$\endgroup$
  • $\begingroup$ 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 at 9:25
  • $\begingroup$ ... there was a question about parsing for differential equation not long ago here pattern-to-match-any-differential-of-a-certain-variable What I am saying, is that if you want to make sure the input is all valid, this is not trivial task. As you have to check for many different combinations and possibilities. $\endgroup$ – Nasser May 14 at 9:28
  • $\begingroup$ @Nasser I am interested in both stages but primarily in the first one, as I expect that the information obtained during the verification stage will make the second stage much easier. $\endgroup$ – jrekier May 14 at 9:33
  • 1
    $\begingroup$ I wrote a small ode solver, just for 1st order ode's and I think 30% of the code was just in parsing. I think WRI internally uses pre-build templates and other internal specialized functions and not how I did his by pattern matching each possibility. In the C days, we used Lex and Yacc to do all the parsing of input. But not now. $\endgroup$ – Nasser May 14 at 9:36
  • $\begingroup$ btw, your last example actually does not work. it gives error from NDSolve. NDSolve::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 at 9:38
9
$\begingroup$

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
   ];
| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for showing these. I did not know about these internal functions. Yes if you could share your long parsing code, that will be good. You do not have to do any customizing for it for this question. I am sure we all will learn from it. Currently for my small ODE solver, I parse things using patterns, but it is a pain and I am always not sure if I missed one case or not and it is only now for just 1 ODE and first order only. $\endgroup$ – Nasser May 14 at 17:32
  • 2
    $\begingroup$ @Nasser Updated. NDSolve`ProcessEquations can be used to parse an IVP. However for a BVP, it will invoke the shooting method to convert it to an IVP, which can be an expensive and problematic process. Then you get all the syntax checking done for you. But it does not seem a good solution for the OP. $\endgroup$ – Michael E2 May 14 at 19:56
  • $\begingroup$ great thanks I will study your code. Very useful. $\endgroup$ – Nasser May 14 at 21:34
  • $\begingroup$ Many thanks. This was exactly the sort of things that I was looking for ! $\endgroup$ – jrekier May 15 at 8:38
7
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you @Nasser, associations look very useful in this context. $\endgroup$ – jrekier May 15 at 8:38
-2
$\begingroup$

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

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

__ WolframAlpha output

This can be driven further:

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

WolframAlpha output

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

WolframAlpha output

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

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ (-1)This question is just not about this, please read it carefully before posting an answer. $\endgroup$ – xzczd Aug 18 at 3:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.