# Improving the speed of evaluating numerical solution to a sytem of PDE

I am running a simulation for a kinetic reaction inside a pellet modeled by a system fo reaction-diffusion equations.

δ = 0.001;
{de, kgco2, kgh2, kgmeoh, kgh2o, kgco} = {1, 1, 1, 1, 1, 1};
r = 8.314;
rho = 8960;


Then, a series of diffusion equations:

e6 = rho*(-r1[x, w] - r2[x, w]) + (de*((2*Derivative[1][psco2][x])/x +
e7 = rho*(-3*r1[x, w] - r2[x, w]) + (de*((2*Derivative[1][psh2][x])/x +
e8 = -(rho*r1[x, w]) + (de*((2*Derivative[1][psmeoh][x])/x +
e9 = rho*(r1[x, w] + r2[x, w]) + (de*((2*Derivative[1][psh2o][x])/x +
e10 = rho*r2[x, w] + (de*((2*Derivative[1][psco][x])/x +


They are accompanied by the following boundary conditions:

ei1 = Derivative[1][psco2][δ] == 0;
ei2 = Derivative[1][psh2][δ] == 0;
ei3 = Derivative[1][psmeoh][δ] == 0;
ei4 = Derivative[1][psh2o][δ] == 0;
ei5 = Derivative[1][psco][δ] == 0;
ei6 = -de*Derivative[1][psco2][1] == kgco2*(psco2[1] - pco2[w]);
ei7 = -de*Derivative[1][psh2][1] == kgh2*(psh2[1] - ph2[w]);
ei8 = -de*Derivative[1][psmeoh][1] == kgmeoh*(psmeoh[1] - pmeoh[w]);
ei9 = -de*Derivative[1][psh2o][1] == kgh2o*(psh2o[1] - ph2o[w]);
ei10 = -de*Derivative[1][psco][1] == kgco*(psco[1] - pco[w]);


With a complicated kinetic function:

r1[x_, w_] := (k5k2k3k4kh2[tp[w]]*psco2[x]*psh2[x]*(1 -
1/kstar[tp[w]]*(psh2o[x]*psmeoh[x])/(psh2[x]^3*psco2[x])))/(1 +
kh2ok8[tp[w]]*psh2o[x]/psh2[x] + rootkh2[tp[w]]*Sqrt[psh2[x]] +
kh2o[tp[w]]*psh2o[x])^3;


As well as:

r2[x_, w_] := (k1[tp[w]]*
psco2[x]*(1 - k3[tp[w]]*(psh2o[x]*psco[x])/(psh2[x]*psco2[x])))/(
1 + kh2ok8[tp[w]]*psh2o[x]/psh2[x] + rootkh2[tp[w]]*Sqrt[psh2[x]] +
kh2o[tp[w]]*psh2o[x]);


The k parameters are defined:

k5k2k3k4kh2[t_] := 1.07*Exp[36696/(r*t)];
rootkh2[t_] := 0.499*Exp[17197/(r*t)];
kh2o[t_] := 6.62*10^-11*Exp[124119/(r*t)];
kh2ok8[t_] := 3453.38;
k1[t_] := 1.22*10^10*Exp[-94765/(r*t)];
kstar[t_] := 10^(3589/t - 11.2);
k3[t_] := 10^(2067/t - 2.0);


With dummy functions:

tp[w_] := 273 + 210;
pco2[w_] := 1;
ph2[w_] := 1;
pmeoh[w_] := 1;
ph2o[w_] := 1;
pco[w_] := 1;


The system is solved through:

sol = ParametricNDSolve[{e6, e7, e8, e9, e10, ei1, ei2, ei3, ei4, ei5,
ei6, ei7, ei8, ei9, ei10}, {psco2, psh2, psmeoh, psh2o,
psco}, {x, δ, 1}, {w}]


Mathematica accepts this system but when I try to obtain a solution:

Plot[psco2[10][x] /. sol, {x, δ, 1}, PlotRange -> All]


The computation takes forever. Is there a way to improve the speed of this computation ?

• w is an internal (dummy) variable in your expressions. It does not appear anywhere in the final equations you pass to ParametricNDSolve. That aside my advice is to first demonstrate that NDSolve works for a particular parameter, before working with ParametricNDSolve. – george2079 Mar 24 '16 at 14:43
• Another suggestion is to try to formulate this as an initial value problem, (all b.c. at delta ). Do that to prove out the consistency of your system, then approach the boundary value problem, which will likely require specifying StartingInitialConditions – george2079 Mar 24 '16 at 15:14
• Try Plot[Evaluate[psco2[10][x] /. sol],....] – user21 Mar 24 '16 at 20:37
• The first steps to improving speed are to eliminate w and replace ParametricNDSolve by NDSolve, replace SetDelayed by Set everywhere, move constants to the beginning of the computation. However, none of these actions seem to have any effect on obtaining an answer in a reasonable time. The root cause, therefore, is the system of ODEs itself. Consider trying the @george2079 's second recommendation. – bbgodfrey Mar 25 '16 at 1:32
• @george2079 I am using ParametricNDSolve because my dummy functions normally depend on w, although I have surpressed the dependence in this example in order to simplify the computation. I'm afraid I don't understand your suggestion about initial value problem, could you elaborate? – Whelp Mar 25 '16 at 9:12

First Attempt

The following is a partial answer. It entails

1. Eliminating w and replaced ParametricNDSolve by NDSolve
2. Replacing SetDelayed by Set everywhere
3. Moving constants to the beginning of the computation
4. Replacing Sqrt[psh2[x]] by Sqrt[Abs[psh2[x]]]
5. Introducing xmax as the upper limit of integration
6. Successfully obtaining an answer for xmax = .1
7. Explicitly using the "Shooting" Method
8. Gradually increasing xmax, using "StartingInitialConditions" obtained from calculations with smaller xmax.

(The first three improve speed modestly but otherwise are not important.) I am stalled at xmax = .75. Here is the code and results.

δ = 0.001; xmax = .75;
{de, kgco2, kgh2, kgmeoh, kgh2o, kgco} = {1, 1, 1, 1, 1, 1};
r = 8.314;
rho = 8960;

tp = 273 + 210;
pco2 = 1;
ph2 = 1;
pmeoh = 1;
ph2o = 1;
pco = 1;

k5k2k3k4kh2[t_] = 1.07*Exp[36696/(r*t)];
rootkh2[t_] = 0.499*Exp[17197/(r*t)];
kh2o[t_] = 6.62*10^-11*Exp[124119/(r*t)];
kh2ok8[t_] = 3453.38;
k1[t_] = 1.22*10^10*Exp[-94765/(r*t)];
kstar[t_] = 10^(3589/t - 11.2);
k3[t_] = 10^(2067/t - 2.0);

e6 = rho*(-r1[x] - r2[x]) + (de*((2*Derivative[1][psco2][x])/x +
e7 = rho*(-3*r1[x] - r2[x]) + (de*((2*Derivative[1][psh2][x])/x +
e8 = -(rho*r1[x]) + (de*((2*Derivative[1][psmeoh][x])/x +
e9 = rho*(r1[x] + r2[x]) + (de*((2*Derivative[1][psh2o][x])/x +
e10 = rho*r2[x] + (de*((2*Derivative[1][psco][x])/x +

ei1 = Derivative[1][psco2][δ] == 0;
ei2 = Derivative[1][psh2][δ] == 0;
ei3 = Derivative[1][psmeoh][δ] == 0;
ei4 = Derivative[1][psh2o][δ] == 0;
ei5 = Derivative[1][psco][δ] == 0;
ei6 = -de*Derivative[1][psco2][xmax] == kgco2*(psco2[xmax] - pco2);
ei7 = -de*Derivative[1][psh2][xmax] == kgh2*(psh2[xmax] - ph2);
ei8 = -de*Derivative[1][psmeoh][xmax] == kgmeoh*(psmeoh[xmax] - pmeoh);
ei9 = -de*Derivative[1][psh2o][xmax] == kgh2o*(psh2o[xmax] - ph2o);
ei10 = -de*Derivative[1][psco][xmax] == kgco*(psco[xmax] - pco);

r1[x_] := (k5k2k3k4kh2[tp]*psco2[x]*psh2[x]*(1 - 1/kstar[tp]*(psh2o[x]*psmeoh[x])/
(psh2[x]^3*psco2[x])))/(1 + kh2ok8[tp]*psh2o[x]/psh2[x] + rootkh2[tp]
*Sqrt[Abs[psh2[x]]] + kh2o[tp]*psh2o[x])^3;
r2[x_] := (k1[tp]*psco2[x]*(1 - k3[tp]*(psh2o[x]*psco[x])/(psh2[x]*psco2[x])))/
(1 + kh2ok8[tp]*psh2o[x]/psh2[x] + rootkh2[tp]*Sqrt[Abs[psh2[x]]] + kh2o[tp]*psh2o[x]);

sol = NDSolve[{e6, e7, e8, e9, e10, ei1, ei2, ei3, ei4, ei5, ei6, ei7,
ei8, ei9, ei10}, {psco2, psh2, psmeoh, psh2o, psco}, {x, δ, xmax},
Method -> {"Shooting", "StartingInitialConditions" -> {ei1, ei2, ei3, ei4, ei5,
psco2[δ] == 1.988032083488816, psh2[δ] == 3.044559836643407,
psmeoh[δ] == 1.31253522848969, psh2o[δ] == 0.4781146784359612,
psco[δ] == 0.13082302200921797}}];

{psco2[x], psh2[x], psmeoh[x], psh2o[x], psco[x]} /. sol /. x -> δ
(* {{2.91349, 4.60682, 1.84667, -0.913494, -0.0668281}} *)

Plot[Evaluate[{psco2[x], psh2[x], psmeoh[x], psh2o[x], psco[x]} /. sol],
{x, δ, xmax}, PlotRange -> All]


Ideally, knowledge of the physical problem being solved would provide good values for "StartingInitialConditions". Even a priori bounds on the "StartingInitialConditions" would be helpful to discourage "Shooting" from sampling the wrong regions of solution space.

Plot[(psh2o[x] - psco[x] + psmeoh[x]) /. sol, {x, δ, xmax}]


reveals that psh2o[x] - psco[x] + psmeoh[x] is equal to 1. with high accuracy. The same is true of psco2[x] + psco[x] - psmeoh[x] and -psh2[x] + -psco[x] + 3 psmeoh[x]. In fact, this is to be expected. A linear combination of the ODEs above yields

(de*((2*Derivative[1][ps][x])/x + Derivative[2][ps][x]))/(r*radius^2*tp) == 0;


where ps represents any of these three linear combinations. When combined with corresponding boundary conditions, it yields ps[x] -> 1. Thus, the ODE system can be simplified significantly by eliminating three of the five dependent variables. With the substitution,

{psh2o[x] -> psco[x] - psmeoh[x], psco2[x] -> -psco[x] + psmeoh[x],
psh2[x] -> -psco[x] + 3 psmeoh[x]}


the functions r1 and r2 become

r1[x_] := (k5k2k3k4kh2[tp] (psco[x] - 3 psmeoh[x]) (psco[x] - psmeoh[x])
(kstar[tp] (psco[x] - 3 psmeoh[x])^3 - psmeoh[x]))/(kstar[tp] (kh2o[tp] psco[x]^2 +
psmeoh[x] (-3 + kh2ok8[tp] + 3 kh2o[tp] psmeoh[x] - 3 Sqrt[Abs[psco[x] - 3 psmeoh[x]]]
rootkh2[tp]) + psco[x] (1 - kh2ok8[tp] - 4 kh2o[tp] psmeoh[x] +
Sqrt[Abs[psco[x] - 3 psmeoh[x]]] rootkh2[tp]))^3)

r2[x_] := (k1[tp] (psco[x] - psmeoh[x]) ((-1 + k3[tp]) psco[x] + 3 psmeoh[x]))/
(kh2o[tp] psco[x]^2 + psmeoh[x] (-3 + kh2ok8[tp] + 3 kh2o[tp] psmeoh[x] -
3 Sqrt[Abs[psco[x] - 3 psmeoh[x]]] rootkh2[tp]) + psco[x] (1 - kh2ok8[tp]
- 4 kh2o[tp] psmeoh[x] + Sqrt[Abs[psco[x] - 3 psmeoh[x]]] rootkh2[tp]))


and the solution is provided by

xmax = 1;
sol1 = NDSolve[{e8, e10, ei3, ei5, ei8, ei10}, {psmeoh, psco}, {x, δ, xmax},
Method -> {"Shooting", "StartingInitialConditions" -> {ei3, ei5,
psmeoh[δ] == 1.2608611256607714, psco[δ] == -0.011072222667300185}}];
{-psco[x] + psmeoh[x], -psco[x] + 3 psmeoh[x], psmeoh[x], psco[x] - psmeoh[x],
psco[x]} /. sol1 /. x -> δ
(* {{1.35351, 4.023, 1.33475, -1.35351, -0.0187626}} *)
Plot[Evaluate[{-psco[x] + psmeoh[x], -psco[x] + 3 psmeoh[x], psmeoh[x],
psco[x] - psmeoh[x], psco[x]} /. sol1], {x, δ, xmax}, PlotRange -> All,
PlotLegends -> {psco2[x], psh2[x], psmeoh[x], psh2o[x], psco[x]}]


That this solution differs from the earlier one (besides the fact that xmax is larger) should not be too surprising. These highly nonlinear equations may have more than one solution. Moreover, I could not extend in the least the earlier solution beyond xmax = 3/4 despite considerable effort. The earlier solution simply may not exist for larger xmax.

• Thanks for taking an interest in this problem. I've been trying to reproduce your results but so far I am unsuccessful. When I naively copy the first code (first attempt) into a blank notebook, I encounter an error message: At x == 0.08293242960370201, step size is effectively zero; \ singularity or stiff system suspected – Whelp Mar 29 '16 at 16:53
• @Whelp Something very strange happening here. I shall post a fix as soon as I can. – bbgodfrey Mar 29 '16 at 18:27
• @Whelp In the meantime, try "Full Solution" To do so, start a fresh session of Mathematica and copy into a new notebook the code from "First Attempt" up to and including ei10. Then, replace xmax=0.75 by xmax = 1 in the first line of the code just copied, Next, copy the last two blocks of code (new r1 and r2 definitions and sol1 = etc.) from "Full Solution". This works for me, and it should for you. Please let me know. Thanks. – bbgodfrey Mar 30 '16 at 0:15
• @Whelp Amazingly, replacing r*radius^2*tp by r*tp*radius^2 solved the problem. Apparently, the calculation in "First Attempt" is on the threshold of instability, and the difference in roundoff is sufficient to change the answer. In any event, I urge you not to waste time on the "First Attempt" solution. "Full Solution" addresses your needs, I believe. – bbgodfrey Mar 30 '16 at 1:36
• Wow, that is very sensitive. Following your steps I've been able to get it to ru, this should be a very strong first step in getting this to work for the full problem. I will now have to introduce temperature dependence as well as axial location in the bed (w). In this kind of numerical simulation with a singularity at x=0, is there a way to avoid setting a δ cutoff? – Whelp Mar 30 '16 at 12:44