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

I have a relatively simple 2nd Order ODE problem I'm trying to solve numerically

$p'' + \frac{1}{r}p' = A(\frac{p(r)}{p(r) + k_{m}})$

where $A = 7.5795*10^9$ and $k_{m} = 1$.

So trying to get numerical solution in Mathematica between the boundaries $r_{o} = 5 \times 10^{-6}$ and $r_{c}$ ($r_{c}$ is at LEAST 0.000102329 and likely 10-30 microns bigger, but I'll use small estimate for this example . We also know $p(r_{o}) = p_{o} = 100$, and $p'(r_{c}) = p(r_{c}) = 0$; but this is seriously upsetting Mathematica; I'm getting the following warnings depending on what versions I run;

s = NDSolve[{D[r*O2'[r], r] - con*r*((O2[r])/(O2[r] + km)) == 0,  O2[rc] == 0, O2[ro] == po}, O2, {r, ro, rc},Method -> {"StiffnessSwitching"}, MaxSteps -> 100000]

Produces a host of errors of the form;

NDSolve::ndsz: At r == 0.00009971646984096604`, step size is effectively zero; >singularity or stiff system suspected

Trying the other condition

s = NDSolve[{D[r*O2'[r], r] - con*r*((O2[r])/(O2[r] + km)) == 0,  O2'[rc] == 0, O2[ro] == po}, O2, {r, ro, rc},Method -> {"StiffnessSwitching"}, MaxSteps -> 100000]

produces the errors

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 >iterations. NDSolve::berr: "There are significant errors {-0.125305,-99.99} in the boundary value >residuals. Returning the best solution found."

The attempted solutions look frankly obscenely wrong, so Im guessing there's a stiffness or singularity probem. I've tried a few tricks to try and circumvent this and I'm getting no joy; is there something I'm missing? So far I've tried increasing the value for rc, and re-writing the equation differently but I'm still having problems. Anyone any idea of why I'm getting these errors and how to resolve them?

share|improve this question
In my understanding, you'd expect a singularity at r = 0, but the function is not evaluated here, and there should be no singular points between ro and rc. Most odd... – DRG Oct 21 '13 at 14:59
up vote 2 down vote accepted

Your equation can be solved by "shooting method" with a slight adjustment of rc:

con = 7.5795 10^9; km = 1; po = 100;
ro = 5 10^-6; rc = 132.329 10^-6;(* Based on your description, I made rc 30×10^-6 bigger. *)
s = NDSolve[{D[r*O2'[r], r] - con*r*((O2[r])/(O2[r] + km)) == 0, O2'[rc] == 0, O2[ro] == po},
            O2, {r, ro, rc}, Method -> {"Shooting", "StartingInitialConditions" -> 
                                                    {O2[rc] == 0, O2'[rc] == 0}}];      
{O2[rc], O2'[rc]} /. s
Plot[O2[t] /. s, {t, ro, rc}, PlotRange -> All]

{{0.167193, 0.}}

enter image description here

O2[rc] is a little different from what you give, I think it doesn't hurt and believe that it can be improved by more adjustment but I don't have time now. (According to my experience, equations which have to turn to shooting method are often sensitive. )

There're many posts about shooting method in this site, you can have a search.

share|improve this answer
That is excellent, thank you - shooting method works fine. As you do point out, it is extremely sensitive; at least this narrows the range I have to look in. Thanks again! – DRG Oct 22 '13 at 9:11

Your Answer


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.