Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
    
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
add comment

1 Answer

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]]

plot

share|improve this answer
add comment

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.