1
$\begingroup$

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 ?

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ 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) + ... $\endgroup$
    – SPPearce
    Commented May 15, 2017 at 12:03
  • $\begingroup$ Just moved the delta away from 0. Would that help ? $\endgroup$
    – Whelp
    Commented May 15, 2017 at 12:34
  • $\begingroup$ It is probably not essential for this set of equations, but it is something to consider if you use higher order differential equations. $\endgroup$
    – SPPearce
    Commented May 15, 2017 at 12:39
  • $\begingroup$ I have seen in the literature that this problem can be solved using the orthogonal collocation method. Is this method implemented in MMA ? $\endgroup$
    – Whelp
    Commented May 15, 2017 at 13:30

1 Answer 1

Reset to default
4
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Whelp
    Commented May 15, 2017 at 12:57
  • $\begingroup$ 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. $\endgroup$
    – Chris K
    Commented May 15, 2017 at 13:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.