current community

your communities

more stack exchange communities

company blog
Stack Exchange Inbox Reputation and Badges
tour help
Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top
up vote 3 down vote favorite

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?

share|improve this question
4  
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}] – Marius Ladegård Meyer 5 hours ago
up vote 5 down vote

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

share|improve this answer
up vote 4 down vote

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!

share|improve this answer
    
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 :) – xzczd 4 hours ago
    
@ xzczd of course, but since this is evident, the OP can figure it out himself and decide, what solution is he interested in. – Alexei Boulbitch 3 hours ago

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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