# 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'. Commented Jan 17, 2016 at 11:01
• The problem is that the set of equations are not a differential equation in $w$. Thus the error message from Mathematica. Commented Jan 17, 2016 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. Commented Jan 17, 2016 at 11:45
• Is ei1 supposed to be bc1, etc.? Commented Jan 17, 2016 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. Commented Jan 17, 2016 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[2][ps][x] + (2/x)*Derivative[1][ps][x] - r[x, w] == 0,
Derivative[1][ps][0] == 0,
Derivative[1][ps][1] == (ps[1] - 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[2]])/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 ? Commented Jan 17, 2016 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. Commented Jan 18, 2016 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/… Commented Jan 18, 2016 at 13:08
• Have you tried my code with one of your real r[x,w] functions? Commented Jan 18, 2016 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. Commented Jan 19, 2016 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'[1] == 0;
bc2[w_] := ps'[2] == ps[2] - 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'[1] == ps[1] - 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. Commented Jan 14, 2016 at 16:41
• This is basically a version of equation 12-11 p. 6 of this document : umich.edu/~essen/html/byconcept/chapter12.pdf Commented Jan 14, 2016 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}$ Commented Jan 15, 2016 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}. >> Commented Jan 15, 2016 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. Commented Jan 17, 2016 at 2:02