# Easy way to extract coefficients of functions in an ODE/PDE? [closed]

Suppose I have an ODE like $$\epsilon y''(x)y(x) + xy'(x) + \epsilon^2 y(x)=0$$ or a PDE like $$\dfrac{{\cal U}}{{\cal T}}\dfrac{\partial u}{\partial t}+\dfrac{{\cal U}^{2}}{{\cal L}}u\dfrac{\partial u}{\partial x}+\dfrac{\epsilon{\cal U}^{2}}{{\cal L}}v\dfrac{\partial v}{\partial y}=0.$$

I had a go at writing some code which specifies what the independent and dependent variables are beforehand, and extracts the part of each of the terms which is not one of those or its derivatives as the coefficient. However, it was a bit haphazard as I often needed to create a table containing all the terms in the equation after rearranging it to have RHS = 0, and then extract the factors from each term in the table that were neither functions nor independent variables.

For the former example, here is some code I have written, with explanations:

1. A function which takes a function and returns all partial derivatives of it, with respect to all its independent variables, up to 4th order (this choice of order is arbitrary)
Derivatives[func_[vars__]] :=
DeleteDuplicates[
Flatten[Table[D[func[vars], {{vars}, k}], {k, 0, 4}]]]

1. Calculation of sets of all derivatives and independent variables in this specific problem.
allDerivatives = Derivatives[y[x]];
allIndepVars = {x};

1. Specification of ODE.
exampleODE = \[Epsilon] y''[x] y[x] + x y'[x] + \[Epsilon] ^2 y[x] == 0;

1. The following: a function which splits apart the factors inside one of the terms in the problem, checks to see if it's an independent variable, dependent variable, or a derivative; a function which returns the product of elements in a list; a function which extracts the non-functional coefficient of a term involving a function or its derivatives; a function which puts the coefficients of all the terms in a list.
NotAVar[expr_] :=
If[Flatten[
Intersection[{expr}, Union[allVars, allDerivatives]]] != {}, 1,
expr]

listProduct[x_List] := Times @@ x;

Coeff[term_] :=
listProduct[Table[NotAVar[term[[k]]], {k, 1, Length[term]}]]

CoeffList[prob_] :=
Table[Coeff[SubtractSides[prob][[1, j]]], {j, 1,
Length[SubtractSides[prob][]]}]


Applying CoeffList[exampleODE] gives the result $$\{\epsilon, 1, \epsilon^2 \}$$, as expected.

Is there a more elegant way of doing it than this?

Update from comment:

For the purposes of my work, polynomial dependencies† can be assumed. I would really like the list of coefficients to be displayed in the same order as their corresponding terms in the equation.

[†I believe that means in the case of an ODE $$F(x,y,y',y'',\dots)=0$$ like the first one above that $$F$$ is a polynomial; likewise for a PDE $$F(t,x,\dots;u,u_t,u_x,\dots)=0$$. — Ed.]

• Please show your expressions as code in addition to the TeX. Don't make us retype your equations in. Also, since you have tried code already, please include that as well. Feb 26, 2021 at 16:41
• @MarcoB Done. Have put a MWE in the question. Feb 27, 2021 at 21:35
• "Applying CoeffList[exampleODE] gives the result {ϵ,1,ϵ^2}, as expected." Your code doesn't give this result, please double check it. Also, what exactly is the expected output? {ϵ,1,ϵ^2}, or {ϵ^2,1,ϵ}? Mar 1, 2021 at 6:15
• Are we to assume that the DE is a polynomial in the dependent variables and their derivatives? Mar 2, 2021 at 19:13
• Does it matter in what order the coefficients are returned? Terms whose coefficient is zero are to be omitted? That is, \[Epsilon] y''[x] + x y'[x] + \[Epsilon] ^2 y[x] == 0 and \[Epsilon] y''[x] y[x] + x y'[x] + \[Epsilon] ^2 y[x] == 0 should return the same thing even though they are different ODEs? Mar 2, 2021 at 19:29

The workflow provided in the OP seems to apply to ordinary differential equations up to order 4. The following workflow is more general, but it will require that parameters and independent variables are expressed as monomials.

First, set up some test ODE's and PDE's:

(*Test ODE's and PDE's*)
odeeqn = ϵ y''[x] y[x] + x y'[x] + ϵ^2 y[x] == 0;
pdeeqn = U/τ D[u[t, x, y], t] +
U^2/L u[t, x, y] D[u[t, x, y], x] + (ϵ U^2)/
L v[t, x, y] D[v[t, x, y], y] == 0;
pdeeqn2 =
U/τ D[u[t, x, y], t] +
U^2/L u[t, x, y] D[u[t, x, y], x] + (ϵ U^2)/
L v[t, x, y] D[v[t, x, y], y] == (Df U)/L^2 D[u[t, x, y], x, x];
pdeeqn3 =
1/L^2 Inactive[
Div][((-κ) .
x], {x}]), {x}] + (ρ Cp)/τ D[Θ[t,
x], t] == 0;
(*Ensure that PDE is of the form LHS\[Equal]0*)
ode = First@SubtractSides@odeeqn;
pde = First@SubtractSides@pdeeqn;
pde2 = First@SubtractSides@pdeeqn2;
pde3 = First@SubtractSides@pdeeqn3;


Set up some helper functions:

(*Extract dependent functions*)
depFns = Union[Cases[#, Derivative[__][f_][v__] -> f[v], All]] &;
(*Extract independent variables*)
indepVbls = Union[Cases[#, Derivative[__][f_][v__] -> v, All]] &;
(*Extract parameters*)
parmFn = Complement[Union[ Cases[#, _Symbol, All]], indepVbls@#] &;
(*Create replacement rules to convert independent variables to 1*)
replFn = Thread[# -> ConstantArray[1, Length[#]]] &;
(*Extract parameter exponent*)
Clear[parmExp]
parmExp[term_][pat_] :=
Quiet@Check[First @ Cases[term, pat^_, All],
Check[First @ Cases[term, pat, All], 1]]
(*Calculate parameter ratios*)
Clear[parmRatio]
parmRatio[op_] :=
Apply[Times, parmExp[#] /@ parmFn@(#)] & /@ List @@ op /.
replFn@indepVbls@op


Now, test the parmRatio function on the test ODE's and PDE's.

Grid[Prepend[{#, parmRatio[#]} & /@ {ode, pde, pde2,
pde3}, {"Operator", "Parameters"}], Frame -> All] The results are consistent with expectations and even appear to work with more complex and inactivated differential forms such as Div and Grad.

Here is an alternate and potentially cleaner way to display the data using datasets.

Clear[opTermParm]
opTermParm[op_] := Module[{terms, parms, ds},
terms = List @@ op;
parms = parmRatio[op];
"Parameters"} -> {##}] & @@@ (Transpose[{List @@ #,
parmRatio[#]}] &@op);
ds = Dataset@ds;
ds[All, <|"Operator" -> op, #|> &]
]
ds = Join @@ (opTermParm[#] & /@ {ode, pde, pde2, pde3});
ds[GroupBy[Key["Operator"] -> KeyDrop["Operator"]]] # Source function extraction

In response to a comment, you may be able to use the following to extract linear source functions that were added to the right-hand side of the ODE. For example:

allFns = Union[Cases[#, _ f_[v_] -> f[v], All]] &;
derivs = (depFns[#]~Join~
Union[Cases[#, Derivative[u__][f_][v__] -> Derivative[u][f][v],
All]]) &;
sourceFns = Union@Complement[allFns[#], derivs[#]] &;


Now, we will test on an example ODE.

odeeqn2 = ϵ y''[x] y[x] + x y'[x] + ϵ^2 y[x] ==
a u[x] + b w[x] + c z[x];
ode2 = First@SubtractSides@odeeqn2;
sourceFns[ode2]
(* {u[x], w[x], z[x]} *)

• Thanks very much for this comprehensive answer. A couple of follow-up questions: - How can I extend the code so that it picks up additional functions which are inside the ODE (e.g. if I added z[x] to the RHS of the first ODE)? - How can I write a command that takes an equation and returns all terms containing a the derivatives of a function as well as those containing the function itself? Mar 11, 2021 at 20:17
• @epsilonD3LT4 I updated the answer addressing your first question (I think). I am not entirely sure what you are asking in your second question. Do you have an example input in an example desired output? Patterns are really not my strong suit so you may be reaching the limits of my capability. Mar 11, 2021 at 21:38
• Consider the equation defined by: pdeeqn = Laplacian[u[x, y], {x, y}] + 2 u[x, y] + u[x, y]^2 + v[x, y] == 0. I've defined a function Terms[eqn_] := (ExpandAll@First@SubtractSides@eqn)[[#]] & /@ Range[Length[ExpandAll@First@SubtractSides@eqn]] which returns all the terms in an ODE/PDE when moved over to one side. I would like to write a function which returns all the terms that are functions of u or functions of the derivatives of u. I have tried Cases[Terms[pdeeqn], (Derivative[n___][u])[v__]*p___ -> (Derivative[n][u])[v]*p, All] but this doesn't work. Mar 11, 2021 at 21:44
• Interestingly, the final command I gave in my last comment does actually work for very complicated partial differential equations (but only extracting terms of functions whose partial derivatives are in the equation) but it doesn't work for the trivial example I give above. However, Total[Cases[Terms[pdeeqn], f_[v__]*p___ -> u[v]*p, All]] gives $2u$, as expected. Is there a way to generalise this idea so that functions of functions can be extracted, so that the $u^2$ term is extracted successfully? Also, is there an equivalent of what I have discussed above, but that Mar 11, 2021 at 23:40