# Method of undetermined coefficients for ODE

### Question

Find the line integral curves of the differential equation $$y' + xy'^2 - y = 0$$.

### Solution:

The line integral curve of the differential equation $$y' + xy'^2 - y = 0$$ is of the form $$y = kx + b$$. Thus, $$y' = k$$. Substituting into the differential equation, we get:

$$k + xk^2 - kx - b = 0$$

From this, we have: $$k=b,$$ $$k=0$$ or $$k=1$$

Therefore, $$k^2=k,$$ $$b=0$$ or $$b=1$$

Thus, the line integral curves of the differential equation $$y' + xy'^2 - y = 0$$ are $$y = 0$$ or $$y = x + 1$$.

This is an exercise from an ODE textbook. The solving process uses the method of undetermined coefficients. I want to solve it automatically using MMA (Mathematica) code. The code is as follows, but I feel it's not very concise. Is there a more elegant and general code?

ClearAll["Global*"];

eq = y'[x] + x*y'[x]^2 - y[x] == 0;

linearForm = k*x + b;

derivativeForm = D[linearForm, x];

reducedEq =
Simplify[eq /. {y[x] -> linearForm, y'[x] -> derivativeForm}];

simplifiedEq = Expand[reducedEq];

constantPart = simplifiedEq /. x -> 0;

xCoefficient = Coefficient[simplifiedEq, x];

kValues = Solve[k - k^2 == 0, k];

associatedSolutions =
Table[{k -> kVal, b -> kVal}, {kVal, kValues[[All, 1]]}]


Additionally, the format of the result obtained from this code is also not good-looking:

{{k->k->0,b->k->0},{k->k->1,b->k->1}}

• "elegant" and "good-looking" are extremely hard to define for code. Can you explain what exactly you would want to optimize in your code? Speed? Length? Readability (although this may also be opinion-based)? Commented Aug 1 at 15:00
• @MarcoB OK. That piece of MMA code is quite similar to doing calculations by hand, and it doesn't really showcase the unique features of MMA. I mean, is there a smarter way to write the code, one that takes better advantage of MMA's automation features, such as built-in functions like FindInstance, ForAll, etc., or perhaps by defining a custom function that is more versatile for solving this type of problem? Commented Aug 1 at 19:04
• This is dAlembert ode. The solution you show is the singular solution. You could used Mathematica DSolve to find it ode=y'[x]+x*y'[x]^2-y[x]==0; sol=DSolve[ode,y[x],x]; sol[[1,2]]/.K[1]->1 Gives y[x] == 1 + x (much easier :) Commented Aug 1 at 23:58
• @Nasser ode=y'[x]+x*y'[x]^2-y[x]==0; sol=DSolve[ode,y[x],x] get: ( $$\text { Solve }\left[\left\{x==\frac{\mathbb{C}_1}{(-\mathbf{1}+\mathrm{K}[\mathbf{1}])^2}+\frac{-\mathrm{K}[\mathbf{1}]+\log [\mathrm{K}[\mathbf{1}]]}{(-\mathbf{1}+\mathrm{K}[\mathbf{1}])^2}, \mathrm{y}[\mathrm{x}]==\mathrm{K}[\mathbf{1}]+\mathrm{xK}[\mathbf{1}]^2\right\},\{\mathrm{y}[\mathrm{x}], \mathrm{K}[\mathbf{1}]\}\right]$$ ) , why K[1]->1? Commented Aug 2 at 0:59
• @lotus2019 Well, Mathematica does not generate singular solution for dAlembert when using IncludeSingularSolutions -> True, (do not know why, it should have, since finding singular solution for dAlembert can be done without finding general solution. So as a workaround I cheated and picked K=1 to make it match what it gave there in sol[[1,2]]. I solved this ode myself and found also $y=0$ to be singular solution as well as $y=1+x$. Commented Aug 2 at 1:15

I do not exactly follow the method your textbook recommend. But this below is how I would solve it and the corresponding Mathematica code below.

## Hand solution

This is standard dAlembert ode, hence it has singular solution(s) and general solution. First will show the hand solution, then the Mathematica code which can be used on any dAlembert type ode to find the singular solutions

Solve $$y^{\prime}+x\left( y^{\prime}\right) ^{2}-y=0$$

Let $$p=y^{\prime}$$

$$p+xp^{2}-y=0$$

Solving for $$y$$ gives $$$$y=xp^{2}+p\tag{1}$$$$

Comparing the above to the standard dAlembert ode form

$$$$y=xf\left( p\right) +g\left( p\right) \tag{2}$$$$

Shows that $$f\left( p\right) =p^{2},g\left( p\right) =p$$. Taking derivative of the above w.r.t $$x$$ gives

\begin{align} y^{\prime} & =f+xf^{\prime}\frac{dp}{dx}+g^{\prime}\frac{dp}{dx}\nonumber\\ p & =f+xf^{\prime}\frac{dp}{dx}+g^{\prime}\frac{dp}{dx}\nonumber\\ p-f & =\left( xf^{\prime}+g^{\prime}\right) \frac{dp}{dx}\tag{3} \end{align}

But $$f^{\prime}=2p,g'=1$$, the above now becomes

$$p-p^{2}=\left( 2xp+1\right) \frac{dp}{dx}$$

Singular solution is when $$\frac{dp}{dx}=0$$ which means $$p-p^{2}=0$$ or $$p\left( 1-p\right) =0$$ or $$p=1,p=0$$. Substituting these into (1) gives the two singular solutions

\begin{align*} y_{1} & =x+1\\ y_{2} & =0 \end{align*}

General solution is when $$\frac{dp}{dx}\neq0$$. This can be solved if needed (not shown).

## Mathematica code

ClearAll["Global*"]
ode = y'[x] + x*y'[x]^2 - y[x] == 0;
ode /. y' -> p;
eq1 = y[x] == First@SolveValues[%, y[x]]


eq3 = Collect[D[eq1, x], p'[x]] /. p'[x] -> 0;
eq3 = eq3 /. y'[x] -> p[x]


pSol = Solve[eq3, p[x]];
singularSols = eq1 /. pSol


• Your approach is more general, more perfect than the solution provided in the textbook, and your code is excellent too. Thank you! Commented Aug 2 at 13:18

SolveAlways[] is the built-in tool for solving for undetermined coefficients. It works well on polynomial-like forms, especially ones in which the form is a linear combination of functions such as in undetermined coefficients in the OP. It can be finicky on other systems (see addendum for an example).

eq = y'[x] + x*y'[x]^2 - y[x] == 0;
form = k*x + b;

coeffs = SolveAlways[eq /. y -> Function @@ {x, form}, x]
sols = y[x] -> form /. coeffs
nontrivialSols = Pick[sols, PossibleZeroQ@Values@sols, False]
(*
{{k -> 0, b -> 0}, {k -> 1, b -> 1}}  <-- coefficients
{y[x] -> 0, y[x] -> 1 + x}            <-- all solutions found
{y[x] -> 1 + x}                       <-- nontrivial solution
*)


This type of example seems to be missing from the docs. The docs have some showing the difference between the generic solutions of SolveAlways and the complete solution of Reduce. None show SolveAlways failing and Reduce succeeding.

SolveAlways[Exp[a x + b] == 2 Exp[2 x], x]
(* < warnings and suggestion to try Reduce omitted >
{{b -> Log[2 E^(2 x - a x)]}}
*)

Reduce[! Reduce[! (Exp[a x + b] == 2 Exp[2 x]), {}, {x}, Reals]]
(*
a == 2 && b == Log[2]
*)


The domain Reals is needed here. Such a restriction cannot be specified for SolveAlways[].

• Your code uses the built-in SolveAlways function really well, making the most of Mathematica's automation features. Plus, you included some special cases of SolveAlways. Great job, and thanks a lot! Commented Aug 3 at 14:23