# Problem solving a second-order PDE

I am looking for numerical solutions for a class of equations of the type :

e = Derivative[2,0][ps][x,w]+(2/x)*Derivative[1,0][ps][x,w]-r[x,w]==0


With r a function that can vary but for the sake of example:

r[x_,w_]:=1-ps[x,w]/2


The boundary conditions are given at x=0 and x=1 by:

bc1=Derivative[1, 0][ps][0, w] == 0
bc2=Derivative[1, 0][ps][1, w] == (ps[1, w] - pc[w])


With pc being taken to be for illustration:

pc[w_] := 5*(1 - 0.01*w)


Mathematica returns an error message upon resolution:

NDSolve[{e, bc1, bc2}, ps, {x, 0, 1}, {w, 0, 10}]


At first I figured I was missing a boundary condition, but adding

bc3=ps[0, 0] == 0


does not seem to help. Any suggestions ?

Edit:

After changing to remove w from NDSolve parameters :

NDSolve[{e, bc1, bc2}, ps, {x, 0, 1}]


The error message comes as :

NDSolve::ndode: Input is not an ordinary differential equation. >>

• You are missing brackets around ps in the equation for e'. – Fabian Jan 17 '16 at 11:01
• The problem is that the set of equations are not a differential equation in $w$. Thus the error message from Mathematica. – Fabian Jan 17 '16 at 11:03
• @Fabian I took your advice and removed the w dependence in NDSolve (ans fixed a typo). I am now getting a different error message, as described in my edit. – Whelp Jan 17 '16 at 11:45
• Is ei1 supposed to be bc1, etc.? – Michael E2 Jan 17 '16 at 13:58
• As I see it, your boundary condition at x=1 depends on the value of ps[1,w]. I don't think NDSolve is going to be happy about that. – Kellen Myers Jan 17 '16 at 14:27

Eliminating w entirely from the ps function (treating it more like a parameter) seems to have some benefit:

r[x_, w_] = 1 - ps[x]/2
pc[w_] = 5*(1 - 0.01*w)
DSolve[
{
Derivative[ps][x] + (2/x)*Derivative[ps][x] - r[x, w] == 0,
Derivative[ps] == 0,
Derivative[ps] == (ps - pc[w])
},
ps,
{x, 0, 1}
]


This returns a function that I've simplified a bit:

f[x_, w_] =
(2. x - 0.0656427 (-60. + 1. w) Sin[x/Sqrt])/x


Unfortunately, DSolve also kicks out the Solve::ifun error, meaning this solution may not be correct.

I've tested the boundary conditions, they seem okay:

Limit[Derivative[1, 0][f][x, w], x -> 0]
(* Output: 0. *)

Chop[Simplify[Derivative[1, 0][f][1, w] - (f[1, w] - pc[w])]]
(* Output: 0. *)


And with some tweaking, I think we can also double check the PDE:

r2[x, w] = 1 - f[x, w]/2;
Chop[FullSimplify[
Derivative[2, 0][f][x, w] + (2/x)*Derivative[1, 0][f][x, w] - r2[x, w]
]]
(* Output: 0. *)


So, with a bit of hacking, I think this settles it?

• So in this case your strategy is to look for an exact solution with w as a parameter, then replace w with its value ? Do you think it will work with more complicated functions for r, with definite dependence on w ? – Whelp Jan 17 '16 at 20:28
• How bad is r[x,w] in general? I mean, if you made it a nightmare, I'm sure DSolve could have issues, but you can change to something like w-ps[x]/2 or even w^2-(w+1)ps[x]/2 and it works (the solution is a bit more of a mess though). In terms of the Mathematica syntax, this strategy seems relatively effective. I can't comment on whether this is theoretically sound in some way -- I'm not especially familiar with this type of PDE. – Kellen Myers Jan 18 '16 at 12:42
• I would say it's pretty bad. In general, r represents the reaction rate and is a complicated function of concentrations and temperatures involving multiple exponentials. A typical expression will look like that :image.slidesharecdn.com/… – Whelp Jan 18 '16 at 13:08
• Have you tried my code with one of your real r[x,w] functions? – Kellen Myers Jan 18 '16 at 13:13
• @Whelp I think treating w as a parameter is a good idea but I don't see why you couldn't do the same with NDSolve. I've added an explanation to my original answer as these comments are character limited. – Philo Jan 19 '16 at 13:53

EDIT:

I think that treating w as a parameter, like @KellenMyers suggested, is a good idea but you should be able to do that with NDSolve as well, not just DSolve, in case your actual problem does not have an analytic solution. For example the following should work (note that I've changed the boundary to be at x=1 and x=2 instead of x=0 and x=1 for reasons that will become apparent)

pc[w_] := 5*(1 - 0.01*w);
r[x_, w_] := 1 - ps[x]/2
e[w_] := ps''[x] + 2/x ps'[x] - r[x, w] == 0;
bc1[w_] := ps' == 0;
bc2[w_] := ps' == ps - pc[w];
sol[w_] := First@NDSolve[{e[w], bc1[w], bc2[w]}, ps, {x, 1, 2}];
f[x_, w_] := ps[x] /. sol[w]


With your example there is an additional problem that one of the boundary conditions is at x=0 and if you try to put that there in the above code you'll get an error of division by 0. I think this happens because the solution is a superposition of spherical bessels whose derivatives diverge at the origin like I pointed out previously. It is only by choosing the coefficients so that the divergent part cancels in the limit $x\rightarrow 0$ that get rid of that. Apparently NDSolve can't handle that. There's probably a fancy and clever way of circumventing this problem but since I'm neither of those the way I would go about it is by shifting your boundary by a tiny amount $x_0$ which is smaller than the accuracy you want to achieve:

pc[w_] := 5*(1 - 0.01*w);
r[x_, w_] := 1 - ps[x]/2
e[w_] := ps''[x] + 2/x ps'[x] - r[x, w] == 0;
bc1[w_, x0_] := ps'[x0] == 0;
bc2[w_] := ps' == ps - pc[w];
sol[w_, x0_] :=
First@NDSolve[{e[w], bc1[w, x0], bc2[w]}, ps, {x, x0, 1+x0}];
f[x_, w_, x0_] := ps[x + x0] /. sol[w, x0]


/EDIT below is my original reply which I later realised was wrong: there is a combination of spherical bessels such that the first derivative doesn't diverge at the origing.

Probably someone else will be able to help you with Mathematica better than I can but, purely on the level of math, it may be worth pointing out that your example is poorly chosen.

Since you don't have any derivatives with respect to w the equation is actually just an inhomogeneous Spherical Bessel equation for the variable x`. The solutions are the 0-order Spherical Bessel functions which diverge at the origin making your boundary conditions impossible.

• Well, the equation is designed to model diffusion/reaction of chemicals inside a catalyst pellet of spherical shape. The boundary conditions describe the absence of particle flow at the center of the pellet (x=0) and at the surface (x=1) the flow is governed by the concentration gradient in the boundary layer. In this, w is a parameter that describe the axial location in a reactor, but the mass transfer I am looking at is only diffusion through the pellet. – Whelp Jan 14 '16 at 16:41
• This is basically a version of equation 12-11 p. 6 of this document : umich.edu/~essen/html/byconcept/chapter12.pdf – Whelp Jan 14 '16 at 16:47
• @Whelp I was wrong, actually only one of the spherical bessels diverges at the origin. The analytic solution to your problem is $$p_s(x,w)=2+\frac{2-p_c(w)}{\frac{\cos(1/\sqrt{2})}{\sqrt{2}}-2\sin(1/\sqrt{2})}\frac{\sin(x/\sqrt{2})}{x}$$ so this should be a well defined problem. Please disregard my previous comment. For some reason LaTeX isn't displayed properly. The last part should be $\frac{\sin(x/\sqrt{2})}{x}$ – Philo Jan 15 '16 at 13:24
• Thanks for the correction. At the moment, I'm still unable to get results. The error mesage reads : NDSolve::femcnmd: The PDE coefficient {{2/x,0}} does not evaluate to a numeric matrix of dimensions {1,2}. >> – Whelp Jan 15 '16 at 15:48
• @Whelp You should include the error message(s) in your question. It would be helpful, especially given all the typos. It's hard to know which problem you were actually facing. – Michael E2 Jan 17 '16 at 2:02