# How does one specify Neumann conditions for NDSolve?

I have a series of functions defined in my notebook, and then want to use this to solve a diffusion-reaction type equation. At the moment, something like this works:

eqn = D[Ef[r, t], t] - Def5*( D[Ef[r, t], r, r] + (2/(r )) D[Ef[r, t], r]) + q[r ] Ef[r, t];

ic = {Ef[r, 0] == kd, Ef[rn, t] == kd, Ef[ro, t] == kd};
s = NDSolve[{eqn == 0}~Join~ic, Ef, {r, rn, ro}, {t, 0, 60}];
Plot3D[Evaluate[Ef[r, t] /. s], {r, rn, ro}, {t, 0, 60},
AxesLabel -> Automatic, MaxRecursion -> 8, PlotPoints -> 32]


q[r] is a predefined function, and rn and ro are boundaries. Now, the above code works, BUT the boundary conditions are wrong and inconsistent; I want to change them to Neumann-type conditions so that $d(Ef)/dr = 0$ at both ro and rn, but I am having trouble implementing this; any ideas? All help appreciated! This is probably very simple but I'm new to Mathematica and sometimes screw up syntax...

DRG

• What kind of trouble are you having exactly? What have you tried? D[Ef[r,t],r]==0 should work Commented Apr 23, 2013 at 13:52
• I've tried something like ic = {Ef[r, 0] == kd, D[Ef[{rn, t}, r]] == 0, D[Ef[{ro, t}, r]] == 0}; but this generates an error "NDSolve::conarg: The arguments should be ordered consistently. "
– DRG
Commented Apr 23, 2013 at 13:53
• If it's as you've written, then the arguments in the boundary conditions are written incorrectly. You have: D[Ef[{rn, t}, r]] whereas you should have D[Ef[rn, t], r] as @Szabolcs has written. Commented Apr 24, 2013 at 2:56
• Very strange; doing that yields the error NDSolve::deqn: Equation or list of equations expected instead of True in the first argument ; really not sure what to make of this...
– DRG
Commented Apr 24, 2013 at 10:42
• ... returning the argument 'true' means it thinks such a statement is self evident, but it surely wouldn't know that unless it evaluates ND solve first?
– DRG
Commented Apr 24, 2013 at 10:53

I got it to work - had to redefine the boundary conditions with Derivative - I'll leave this up here if it's useful to others - working code block is

ic = {Ef[r, 0] == kd, Derivative[1, 0][Ef][rn, t] == 0,
Derivative[1, 0][Ef][ro, t] == 0};

• Thanks for posting the solution, it's good for the community to do that. Note that D[Ef[rn,t], rn] also returns Derivative[1, 0][Ef][rn, t]. I prefer the former because it explicitly mentions the variable according to which the differentiation is done (Derivative doens't). Commented Apr 24, 2013 at 12:41

The solution is either to write the derivative as suggested which first computes the derivative by the first variable and then plugs in the constant value instead of the variable, i.e
Derivative[1,0][ef][rn,t]
or you can write the equivalent by first computing the derivative and then setting the boundary value as
D[ef[r, t], r] /. r -> rn
That would be explicit in the variable name (may require putting in brackets if part of a complex statement as otherwize everything after the -> is used)

What Szabolcs says is not true as rn is a number so in this case
D[ef[rn,t],rn] evaluates to something like D[ef[0,t],0] which doesn't make sense
D[ef[rn,t],r] Doesn't work as well since ef is not a function of r so the derivative is 0