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 solve a differential equation numerically. So I have

zSolutionRule = NDSolve[
   {
    z''[x]*z[x] - z'[x]^2 - z[x]^2 - p0[x]*z[x]^3 == 0,
    p0'[x] == 0,
    z'[0]*E == z'[1],
    z[0] == 1,
    z[1] == E
   },
   {z, p0}, {x, 0, 1}
  ];

When I run this, though, I get an "infinite expression $\frac{1}{0}$ encountered" error, but I want to get the solution to the differential equation. How do I get NDSolve to give me the solution like it does with other differential equations? (I am using Mathematica 9.0.0.0.)

Background

At first I solved a related differential equation

ySolutionRule = NDSolve[
   {y''[x] - 1 - p0[x]*Exp[y[x]] == 0,
    p0'[x] == 0,
    y'[0] == y'[1],
    y[0] == 0,
    y[1] == 1
    },
   {y, p0}, {x, 0, 1}
   ];
ySolution = y /. ySolutionRule[[1, 1]];
chargeNormalization = (p0 /. ySolutionRule[[1, 2]])[.5];

and this worked just fine, but I figured I could try to change variables to see if eliminating the exponential would change how good the solution is if I change boundary conditions when y[1] is big.

So I rewrote the equation in terms of $z=e^y$. I can check that the exponential of ySolution satisfies the transformed differential equation and boundary conditions:

zSol[x_] := Exp[ySolution[x]];
zSol[1]/zSol[0] - E
zSol[0] - 1
zSol[1] - E

and

Plot[zSol''[x]*zSol[x] - zSol'[x]^2 - zSol[x]^2 -chargeNormalization*zSol[x]^3, {x, 0, 1}]

It does, up to some accuracy. So then I tell it to find my solution but instead, it gives me an infinite expression error.

share|improve this question
1  
your p0 is a constant ... –  belisarius Apr 29 '13 at 17:46
    
yes p0 is a constant because I have the equation p0'[x]==0. However, I don't know what the value of the constant is before I call NDSolve; NDSolve is supposed to figure that out for me. I think you are with me in wishing there were a better way of telling NDSolve that p0 was just a constant, but I couldn't find one. I think there might not be one because the pattern I used is given in the documentation: [reference.wolfram.com/mathematica/tutorial/… –  NowIGetToLearnWhatAHeadIs Apr 29 '13 at 23:07
add comment

1 Answer

up vote 4 down vote accepted

By default, NDSolve uses a "shooting" method to satisfy the boundary conditions: it picks an initial condition and then evolves the equation to see what boundary values are produced by the initial conditions. It then rejiggers these initial conditions to produce the boundary values specified by the boundary conditions. Here is the Mathematica documentation for the "Shooting" Method. It turns out that, for my differential equation, NDSolve's initial guess for initial condition is a bad guess. The documentation says

The shooting method by default starts with zero initial conditions so that if there is a zero solution, it will be returned.

This should sound bad because $z=e^y$ so z[0]==0 does not seem like a good guess. So just tell NDSolve what initial guess to make:

zSolutionRule = NDSolve[
   {
    z''[x]*z[x] - z'[x]^2 - z[x]^2 - p0[x]*z[x]^3 == 0,
    p0'[x] == 0,
    z'[0]*E == z'[1],
    z[0] == 1,
    z[1] == E
    },
   {z, p0}, {x, 0, 1},
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {z[0] == 1, z'[0] == 1, 
   p0[0] == -.1}}
   ];
zSolution = z /. zSolutionRule[[1, 1]];
chargeNormalization = (p0 /. zSolutionRule[[1, 2]])[.5];

The above code does not give an error. We can verify that it gives a solution satisfying the boundary conditions:

zSolution[1]/zSolution[0] - E
zSolution[0] - 1
zSolution[1] - E

and

Plot[zSolution''[x]*zSolution[x] - zSolution'[x]^2 - zSolution[x]^2 - 
  chargeNormalization*zSolution[x]^3, {x, 0, 1}, PlotRange -> All]
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.