# How to determine if RHS of first order ODE is separable or linear or neither? Parsing problem

## Introduction

This is really a parsing question. But the context is that given the right side of first order ode, I want to automatically determine the class of the RHS. It can be either separable, or linear. If it is none of these, then I can do something else.

For learning, I am writing a step-by-step ODE solver in Mathematica.

Currently the correct solver is explicitly called based on type of RHS, by first visually looking at it and seeing which form the RHS is.

So I have solver for linear RHS, one for separable RHS, one for exact RHS, etc...

But I want to automate this part, just like one calls DSolve and have the program figure which solver to call based on the type of the RHS.

But the problem it is very hard to parse the RHS and determine the form. I do not know how DSolve does all of this.

## description of problem

Given a first order ODE in $x$ and $y(x)$, the general form is $y'(x)=f(x,y(x))$

So need a function, which takes in $f(x,y(x))$ and return if it is linear, or separable or neither and also return back these parts.

$f(x,y(x))$ is separable if it can be written as product of two functions, where one function can depend on $x$ but not on $y(x)$ and the other function does not depend explicitly on $x$ but can depend on $y(x)$.

As in $f(x,y(x)) = p(x) q(y(x))$ where $p(x)$ can depend on $x$ but not on $y(x)$ and $q$ can depend on $y(x)$ but not on $x$. Examples of separable are

$f=x^2+x^2 y(x) = x^2 (1+y(x))$. Here $p=x^2$ and $q=1+y(x)$

$f=1+\cos(x)$ Here $p=1+\cos(x)$ and $q=1$

$f=y(x)$ Here $p=1$ and $q=y(x)$

There are some which are tricky, such as

$f=1−x+y^2(x)−x y^2(x)$ this is separable since it can be written as $(1−x)(1+y^2(x))$ and in this case $p=(1−x)$ and $q=1+y^2(x)$

Also $f=1$ is separable. Here $p(x)=1$ and $q(y(x))=1$.

$f(x,y(x))$ is linear if it can be written as sum of two functions as follows $f(x,y(x))=p(x)+y(x) q(x)$ where $p(x)$ and $q(x)$ here can depend only on $x$.

Examples of a linear $f$ is $f=x^2+x^3 y(x)$, here $p=x^2$ and $q=x^3$. This can't be written as separable. Also $f=1+ x y$ is linear where $p=1$ and $q=x$. This can't be written as separable.

But $f$ sometimes can be both linear and separable at same time. In this case, I want to classify it as separable. For example $f=x^2+x^2 y(x)$ can be viewed as $x^2(1+y(x))$ which is separable, or as linear where $f=p(x)+y(x)q(x)$ with $p=q=x^2$

Some RHS is not separable nor linear. These I will just flag as neither and worry about them later. For example $f=\frac{1}{x+y}$ as this can't be written in the form $p(x)+y(x)q(x)$ and it also can't be written as $f=p(x) q(y(x))$

## Input for testing

Here are test inputs, and expected output to make it easier to answer this question.

input={
x^2+x^2 y[x],
2 x+1,
Cos[2 x],
y[x],
(x-1)/y[x],
Log[1+y[x]^2],
1-x+y[x]^2-x y[x]^2,
x^2+x^3 y[x],
1/(x+y[x]),
(1+Sqrt[x])/(1+Sqrt[y[x]]),
(1-x^2+y[x]^2-x^2 y[x]^2)/x^2,
2 x y[x]^2+3 x^2 y[x]^2,
3 Exp[2 x]+2 y[x],
(5 Sqrt[x]-y)/ x,
2 x y[x]+3 x^2 Exp[x^2],
Sqrt[1+x+y[x]^2],
1,
1+x y[x]
}


And expected output for each one these inputs is

out={
{x^2,1+y[x],"separable"},
{2x+1,1,"separable"},
{1,y[x],"separable"},
{Cos[2 x],1,"separable"},
{(x-1),1/y[x],"separable"},
{1,Log[1+y[x]^2],"separable"},
{(1-x),1+y[x]^2,"separable"},
{x^2,x^3,"linear"},
{"","","neither"},
{(1+Sqrt[x]),1/(1+Sqrt[y[x]]),"separable"},
{(1-x^2)/x^2,1+y[x]^2,"separable"},
{2 x+3 x^2,y[x]^2,"separable"},
{3 Exp[2 x],2,"linear"},
{5 Sqrt[x]/ x,-1/x,"linear"},
{3 x^2 Exp[x^2],2 x,"linear"},
{"","","neither"},
{1,1,"separable"},
{1,x,"linear"}
}


This is summarized in

title={"f(x,y(x)","p","q","class"}
Grid[Join[{title},Table[{input[[i]],out[[i,1]],out[[i,2]],out[[i,3]]},
{i,Length@out}]],Frame->All]


For example, let the parsing function be called parse then the call using the first example is

parse[ x^2+x^2 y[x], x, y[x]]


And this should return

 {x^2, 1 + y[x], "separable"}


as in

parse[f_, x_, y_] := Module[{p, q, class},
If[FreeQ[f, y[x]],
p = f; q = 1; class = "seperable";
Return[{p, q, class}, Module]
];

{p, q} = f /. Times[p_, q_] :> {p, q};
(*etc....not working*)
]


Basic parsing methods do not really work. For example doing something like

expr=x^2+x^2 y[x]
{p,q}=Simplify[expr]/.Times[p_,q_]:>{p,q}


Works, but not here

expr=1-t+y[t]^2-t y[t]^2
expr/.Times[p_,q_]:>{p,q}


Using Simplify sometimes help and sometimes now. So this approach is not robust to do what I want.

How does DSolve manages to parse such complicated input and it figures which part of the expression is which? I am sure it does something much more advanced than basic pattern matching.

What is the general approach to write such a parser? Can just using pattern matching function work? Or it needs special parsing methods? May be a hint towards a general approach to solve this will be enough. Since trying to do this case by case is not getting to anywhere.

• for linear part of the question, perhaps InternalLinearQ[#, y[x]] & /@ input? – kglr Dec 21 '17 at 9:13
• For the "separable" part, you could use my answer (using separate[expr, {x, y[x]}] or @LeonidShifrin's answer (using getGX). – Carl Woll Dec 21 '17 at 16:07

Thanks for the comments above, I was able to make the first basic working example. This checks for separable parts with the help of code by Leonid Shifrin posted at how-can-i-separate-a-separable-function and the linear part thanks to hint by kglr above using InternalLinearQ and also thanks to CarlWoll for his pattern matching function getPatterns

input={x^2+x^2 y[x],2 x+1,Cos[2 x],y[x],(x-1)/y[x],
Log[1+y[x]^2],1-x+y[x]^2-x y[x]^2,x^2+x^3 y[x],1/(x+y[x]),
(1+Sqrt[x])/(1+Sqrt[y[x]]),(1-x^2+y[x]^2-x^2 y[x]^2)/x^2,
2 x y[x]^2+3 x^2 y[x]^2,3 Exp[2 x]+2 y[x],(5 Sqrt[x]-y[x])/x,
2 x y[x]+3 x^2 Exp[x^2],Sqrt[1+x+y[x]^2],1,1+x y[x]};

out=(dsolveStep[y'[x]==#,y[x],x] &)/@input;

title={"f(x,y(x)","p","q","class"};
Grid[Join[{title},Table[{input[[i]],out[[i,1]],out[[i,2]],out[[i,3]]},
{i,Length@out}]],Frame->All]


Here is the code. This is beta version.

(*This function thanks to  Carl Woll, see
https://mathematica.stackexchange.com/questions/151850/using-cases-\
and-when-to-make-input-a-list-or-not
*)
getPatterns[expr_, pat_] :=
Last@Reap[expr /. a : pat :> Sow[a], _, Sequence @@ #2 &];

(*Thanks to Leonid Shifrin
https://mathematica.stackexchange.com/questions/8142/how-can-i-\
separate-a-separable-function**)
ClearAll[getGX];
getGX[expr_, xvar_, yvar_] :=
With[{dlogg = D[Log[expr], xvar] // FullSimplify},
Exp[Integrate[dlogg, xvar]] /; FreeQ[dlogg, yvar]];

Clear[getHY];
getHY[expr_, xvar_, yvar_] :=
FullSimplify[(#/getGX[#, xvar, yvar]) &[expr]]

(*-------- special case , RHS is constant -----*)
processRHSConstant[ode_, y_[x_], x_] := Module[{},
Print["In processRHSConstant"]
];

(*------- error checking -------------*)
errorChecking[odeInput_, y_[arg_], x_] :=
Module[{ode, lhs, rhs, order, tmp, maxOrder, n, z, independentVariables},

(*check if y[x] is function of x*)
If[Not[SameQ[arg, x]],
Return[Row[{"Argument of dependent variable ", y, " must be ", x}], Module]];

(*check if ODE is equation*)
Return[Row[{"Expected equation as input but found ", odeInput}], Module]];

(*check order of ODE. We only do first order at this time *)
tmp = getPatterns[odeInput, Derivative[n_][y][x]];
order = Cases[tmp, Derivative[n_][y][x] :> n];
If[order === {},
Return[Row[{"No ", y'[x], " found in the ODE ", odeInput}], Module]];

maxOrder = Max[order];
If[maxOrder != 1,
Return[Row[{"Only first order is supported at this time. Found ODE of order ", maxOrder}], Module]];

(*check that dependent variables all depend on x*)
tmp = getPatterns[odeInput, y[z_]];
independentVariables = Cases[tmp, y[z_] :> z];
If[Length[independentVariables] > 0 &&
Not[AllTrue[independentVariables, SameQ[#, x] &]],
Return[Row[{"dependent variable ", y, " must be function of ", x, " only"}], Module]
];

{}
];

(*------- main function -------*)
dsolveStep[odeInput_, y_[arg_], x_] :=
Module[{status, ode, lhs, rhs, tmp, p, q, class, any},

status  = errorChecking[odeInput, y[arg], x];
If[Not[status === {}], Return[status, Module]];

(*convert to standard for y'(x) = f(x,y(x)) *)
lhs = odeInput /. lhs_ == rhs_ :> lhs;
rhs = odeInput /. lhs_ == rhs_ :> rhs;
rhs = -(lhs /. (y'[x] + any_.) :> any) + rhs;
ode = y'[x] == rhs;

(*special cases*)
If[FreeQ[rhs, x], processRHSConstant[ode, y[x], x]];

(*check if separable*)
tmp = getGX[(rhs /. y[x] -> y), x, y];

p = tmp /. y -> y[x];
q = getHY[(rhs /. y[x] -> y), x, y];
q = q /. y -> y[x];
class = "separable"
,(*not separable check linear*)
If[InternalLinearQ[rhs, y[x]],
q = getPatterns[rhs, any_. *y[x]];
q = First@Cases[q, any_. *y[x] :> any];
p = rhs - (q*y[x]);
class = "linear"
,
p = {};
q = {};
class = "neither"
]
];
{p, q, class}
]
`