4
$\begingroup$

For an ODE like this:$(1-y)y'+y^2=0$ with the initial condition $y(1)=1$, how to solve it numerically? I know this equation can be solved analytically by DSolve. In fact, my equation is more complicated than this, I have to solve it numerically. Using NDSolve directly,

NDSolve[{(1 - y[x])*y'[x] + y[x]^2 == 0, y[1] == 1}, y, {x, 1, 5}]

it will display error messages:

Power::infy: Infinite expression 1/0. encountered.
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 1.`.

I guess this problem happens because the initial condition just makes the coefficient of y'[x] be zero. So my question is how to overcome this problem?

$\endgroup$
  • 4
    $\begingroup$ The derivative diverges at the point $x=1$, and it's hard to deal with divergencies numerically. Perturbing the initial condition slightly removes the problem, e.g. NDSolve[{(1 - y[x])*y'[x] + y[x]^2 == 0, y[1] == 1.01}, y, {x, 1, 5}] $\endgroup$ – Marius Ladegård Meyer Nov 22 '16 at 9:57
5
$\begingroup$

One can here introduce another dependent variable: z[x]->y[x] - y[x]^2/2and express your equation in terms of this variable:

    ss = NDSolve[{z'[x] + (1 - Sqrt[1 - 2 z[x]])^2 == 0, z[1] == 1/2}, 
       z[x], {x, 1, 10}][[1, 1]]
(*  z[x] -> InterpolatingFunction[{{1., 10.}}, <>][x]  *)

which is nicely solved:

Plot[1 - Sqrt[1 - 2 z[x]] /. ss, {x, 1, 10}, 
 AxesLabel -> {Style["x", 18, Italic], Style["y", 18, Italic]}]

yielding

enter image description here

Have fun!

$\endgroup$
  • $\begingroup$ This transformation actually will lead to 2 equations: ((1 - y[x])*y'[x] + y[x]^2 == 0 /. Solve[z[x] == y[x] - y[x]^2/2, y[x]] /. (y[x] -> a_) :> (y -> (Function[x, #] &@a))), which correspond to 2 solutions of the original equation :) $\endgroup$ – xzczd Nov 22 '16 at 11:12
  • $\begingroup$ @ xzczd of course, but since this is evident, the OP can figure it out himself and decide, what solution is he interested in. $\endgroup$ – Alexei Boulbitch Nov 22 '16 at 12:26
6
$\begingroup$

Your specific example is actually special, it has 2 solutions, and DSolve can only find one of them. To find both of the solutions, we can modify the equation from a equation of $y(x)$ to a equation of $x(y)$:

$$\frac{1-y}{x'(y)}+y^2=0$$

Then DSolve and NDSolve can both handle the equation without difficulty:

asolinverse = x /. First@DSolve[{(1 - y)*1/x'[y] + y^2 == 0, x[1] == 1}, x, y]
(* Function[{y}, (1 + y Log[y])/y] *)
nsolinverse = 
 x /. First@NDSolve[{(1 - y)*1/x'[y] + y^2 == 0, x[1] == 1}, x, {y, 10^-3, 100}]

ParametricPlot[{nsolinverse[y], y}, {y, 10^-3, 10}, AspectRatio -> 1/GoldenRatio]

Mathematica graphics

$\endgroup$
  • $\begingroup$ your method is also very good. But I think Boulbitch's method is more universal. For more complicated ODEs, it is not easy to modify the equation from an equation of y(x) to an equation of x(y) . $\endgroup$ – Mark_Phys Nov 23 '16 at 10:46

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.