Solving a diffusion-reaction ODE with NDSolve

I'm working on solving an equation to estimate production in a spherical catalyst pellet. A working example is given below:

delta = 0.01;
k=1;
t0 = 260 + 273;
rt = 8.314;
presList = {p[x], p2[x], p3[x], p4[x]};
presBCList = 80*{9/100, 75/100, 0.17/100, 0.07/100};
{a1, a2, a3, a4, a5} = {4224.40, 0.00001446, 3453.38, 24.18, 96.629};

p2[x_] = presBCList[[2]] - 3 (presBCList[[1]] - p[x]);
p3[x_] = presBCList[[3]] + 1 (presBCList[[1]] - p[x]);
p4[x_] = presBCList[[4]] + 1 (presBCList[[1]] - p[x]);


The reaction speed is given by

r[x_] = k ((a1*p[x]*
p2[x]*(1 - 1/a2*(p3[x]*p4[x])/(p2[x]^3*p[x])))/(1 +
a3*p4[x]/p2[x] + a4*Sqrt[p2[x]] + a5*p4[x])^3);


Next setting up NDSolve with no flux at pellet center BC and a given flux at the surface (x=1):

 diffEqList = {(Derivative[2][p][x] + 2/x Derivative[1][p][x]) -
r[x] == 0};
bcList = {Derivative[1][p][delta] == 0,
Derivative[1][p][1] == -(p[1] - presBCList[[1]])};
sol = NDSolve[Join[diffEqList, bcList], p[x], {x, delta, 1}];


The equation is immediatly solved when k=1 to k=11. However, when k>= 11.95, NDSolve fails.

Power::infy: Infinite expression 1/0. encountered.


Yet looking at the graphs, I can not see what quantity may be going to 0. As far as I can tell, the solutions should not be qualitatively different when k is increased. Changing the value of delta to increase it seems to help a little, but not much. Can this equation be solved with a better method ?

• As well as Chris's answer, are you expanding the spherical coordinates around x=0 to get the boundary condition there, or just moving it delta away from 0? i.e. r'(delta) = r'(0) + delta r''(0) + ... May 15, 2017 at 12:03
• Just moved the delta away from 0. Would that help ? May 15, 2017 at 12:34
• It is probably not essential for this set of equations, but it is something to consider if you use higher order differential equations. May 15, 2017 at 12:39
• I have seen in the literature that this problem can be solved using the orthogonal collocation method. Is this method implemented in MMA ? May 15, 2017 at 13:30

NDSolve is using the shooting method for this boundary value problem. To help get it on the right track, you can manually give initial guesses. As noted by @Szabolcs in this answer "this is documented in the Advanced Numerical Differential Equation Solving tutorial".

This seems to help with your problem:

k = 15;
r[x_] = k ((a1*p[x]*p2[x]*(1 - 1/a2*(p3[x]*p4[x])/(p2[x]^3*p[x])))/(1 + a3*p4[x]/p2[x] + a4*Sqrt[p2[x]] + a5*p4[x])^3);

diffEqList = {(Derivative[2][p][x] + 2/x Derivative[1][p][x]) - r[x] == 0};
bcList = {Derivative[1][p][delta] == 0, Derivative[1][p][1] == -(p[1] - presBCList[[1]])};
sol = NDSolve[Join[diffEqList, bcList], p, {x, delta, 1},
Method -> {"Shooting", "StartingInitialConditions" -> {p[delta] == 6, p'[delta] == 0}}];


works fine for me.

• This seems to help a lot indeed. A basic question but how do you get values or plot p'[x] ? I'm not able to find the right combination of commands to do it. Plot[Evaluate[p'[x] /. sol], {x, delta, 1}] does not seem to work. May 15, 2017 at 12:57
• Better to make the unknown in NDSolve p, not p[x]. I'll edit my answer to reflect this. Then Plot[p'[x] /. sol[[1]], {x, delta, 1}] works. May 15, 2017 at 13:40