Singularity / stiffness / precision errors with NDSolve

Question

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?

Answer 1

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.}}

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.