# Differential equations: Solving a singular/nonlinear ODE using NDSolve

I am trying to numerically solve a rather horrible looking singular/nonlinear ode: A*(y''[x]/y[x])'' +y[x]^2==1 on the interval [0,1], where A is a small parameter (~10^(-4)). The boundary conditions are :

y[0]=y[1]==1
y'[0]==0
y''[1]==0.


I expand the term in parentheses and use NDSolve:

s = NDSolve[{10^(-4)*(y''''[x]/y[x] - 2*y'''[x]*y'[x]/(y[x]*y[x]) - (y''[x]/y[x])*
(y''[x]/y[x]) + 2*(y''[x]*y'[x]*y'[x])/(y[x]*y[x]*y[x])) + y[x]*y[x] == 1,
y[0] == 0, y[1] == 0, y'[0] == 0, y''[1] == 0}, y, {x, 0, 1}]


And I get the following errors:

Power::infy: Infinite expression 1/0.^2 encountered. >>
∞::indet: Indeterminate expression 0. ComplexInfinity encountered. >>
Power::infy: Infinite expression 1/0.^2 encountered. >>
∞::indet: Indeterminate expression 0. ComplexInfinity encountered. >>
Power::infy: Infinite expression 1/0.^3 encountered. >>
General::stop: Further output of Power::infy will be suppressed during
this calculation. >>
∞::indet: Indeterminate expression 0. ComplexInfinity encountered. >>
General::stop: Further output of ∞::indet will be suppressed during
this calculation. >>
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0..


I realize the solution to this function is very ill behaved near $x=0$. I was wondering whether there is a way to over come this problem? I know that, in the interior, y[x] = -1.

Thanks a lot!

-
According to Maple, the system appears to be singular at the left point. It actually solved it analytical, with the trivial solution $y=0$. screen shot !Mathematica graphics (so basically your system is not well formed to solve) it suggested a mid-point numerical method. But I do not know now how to apply this in Mathematica NDSolve and it is there... –  Nasser Oct 22 '13 at 0:06
@Nasser I'd say that y[x_]=1 is a better solution. –  b.gatessucks Oct 22 '13 at 15:36

First off, I wanted to see if I could get the solution y[x_] = -1 as claimed. It seems that I can't.

eq = bigA D[(y''[x]/y[x]), {x, 2}] + y[x]^2 - 1 ;
cond = {y[0] == 1, y[1] == 1, y'[0] == 0, y''[1] == 0} ;


Our test function is a polynomial :

test = a0 + a1 # + a2 #^2 + a3 #^3 + a4 #^4 &;


The initial/boundary conditions give :

r = Reduce[cond /. y -> test] // ToRules
(* {a3 -> -((5 a4)/2), a2 -> (3 a4)/2, a1 -> 0, a0 -> 1} *)


We can get the last coefficient by imposing that the equation be satisfied for every x :

last = SolveAlways[((eq /. y -> test) /. r ) == 0, x]
(* {{a4 -> 0}} *)


And finally :

test[x] /. r /. First@last // Simplify
(* 1 *)


Next we can use this insight and make a simple transformation; it turns out that NDSolve will return an answer without complaining (but no exciting solution either).

sol = f /.
First@NDSolve[{(eq /. y -> (f[#] + 1 &) ) == 0 /. bigA -> 1/1000,
cond /. y -> (f[#] + 1 &) }, f, {x, 0, 1}] ;

Plot[sol[x], {x, 0, 1}, PlotStyle -> Thickness[0.02]]
`

-