5
$\begingroup$

I have the following code

ep = 80;
a = 10;
cb[1] = 0.001;
cb[2] = 0.001;
zf = 10;
z[1] = 1;
z[2] = -1;
cp[i_] = cb[i] Exp[-z[i] f[x]];

p = 
  NDSolve[
    {f''[x] == D[HeavisideTheta[x - 10], x] zf/(ep a^2) - 
                 1 Sum[z[i] cp[i], {i, 1, 2}], 
     f[10] == zf/(ep a^2), f'[20] == 0},
     f, {x, 10, 100 }, 
     MaxSteps -> Infinity, 
     InterpolationOrder -> All, 
     Method -> 
       {"Shooting", 
          "StartingInitialConidtions" -> {f[10] == zf/(ep a^2), f'[20] == 0}}]

When I execute the code, I get the following error message:

NDSolve::bvdisc: NDSolve is not currently able to solve boundary value problems with discrete variables. >>

Could you please tell me where I have made a mistake?

$\endgroup$
  • 2
    $\begingroup$ first thing you should do is verify you can solve an lnitial value problem (specify f'[10]) before attempting to use the shooting method. also for the shoting method the StartingInitialConditions are specified at the initial point (10) $\endgroup$ – george2079 Dec 31 '16 at 14:34
  • 2
    $\begingroup$ The problem lies with the delta function, D[HeavisideTheta[x - 10], x] , but why NDSolve chooses to respond with this particular error message is not obvious to me. $\endgroup$ – bbgodfrey Dec 31 '16 at 14:40
  • $\begingroup$ Possible duplicate: (55986). Also related: (133493) $\endgroup$ – Michael E2 Jan 1 '17 at 4:20
5
$\begingroup$

Actually, simply

p = First@NDSolve[{f''[x] == - Sum[z[i] cp[i], {i, 1, 2}], 
    f[10] == zf/(ep a^2), f'[20] == 0}, f, {x, 10, 100 }];
Plot[(f /. p)[x], {x, 10, 100 }]

enter image description here

gives the same answer. The term D[HeavisideTheta[x - 10], x] zf/(ep a^2) has no effect, because it is at the boundary. It is strange that NDSolve gives an error instead of just ignoring the term.

$\endgroup$
  • $\begingroup$ Thanks alot, yes , the problem for me is why mathematica does not understand this. $\endgroup$ – amin bk Dec 31 '16 at 19:48
  • 1
    $\begingroup$ @aminbk it is difficult to say why Mathematica behaves as it does. I recommend that you focus on working around the problem, as I did in this particular case. Do you have a more general problem that you are trying to solve? $\endgroup$ – bbgodfrey Dec 31 '16 at 19:56
4
$\begingroup$

DiracDelta[x] can be replaced by:E^(-(Abs[x]/epsilon))/(2 epsilon) where epsilon is a sufficiently small real number,but not strictly zero.

p = With[{eps = 1/10^9}, 
NDSolve[{f''[x] == E^(-(Abs[-10 + x]/eps))/(2 eps)*zf/(ep a^2) - 
Sum[z[i] cp[i], {i, 1, 2}], f[10] == zf/(ep a^2), f'[20] == 0},f, {x, 10, 100}]];

Plot[f[x] /. p, {x, 10, 100}, AxesOrigin -> {10, 0}]

enter image description here

$\endgroup$
  • $\begingroup$ Why the solution curve starts in the second quadrant? $\endgroup$ – zhk Dec 31 '16 at 16:07
  • 1
    $\begingroup$ @MMM. Problem is a Plot gives strange output.I'm edited answer. $\endgroup$ – Mariusz Iwaniuk Dec 31 '16 at 16:21
4
$\begingroup$

This is more like a comment.

You have a discontinuous ode, so to turn off the discontinuity during NDsolve's processing of your ode with DiscontinuityProcessing,

NDSolve[{f''[x]==DiracDelta[-10 + x] zf/(ep a^2) - 1 Sum[z[i] cp[i], {i, 1, 2}],
f[10] == zf/(ep a^2), f'[20] == 0}, f, {x, 10, 100},
Method -> {"DiscontinuityProcessing" -> False}]

which generate this error,

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 10.`.

So I added a submethod Method -> "ExplicitEuler"

p = NDSolve[{f''[x] == 
  DiracDelta[-10 + x] zf/(ep a^2) - 1 Sum[z[i] cp[i], {i, 1, 2}], 
  f[10] == zf/(ep a^2), f'[20] == 0}, f, {x, 10, 100}, 
  Method -> {"FixedStep", Method -> "ExplicitEuler", 
  "DiscontinuityProcessing" -> False}]

which produced a solution but with a warning,

NDSolve::nlnum: The function value {-0.0000234567,2.5*10^-6+0.00125 DiracDelta[0.]} is not a list of numbers with dimensions {2} at {x,f[x],(f^[Prime])[x]} = {10.,0.00125,-0.0000234567}.

Plot[f[x] /. p, {x, 10, 100}]

enter image description here

Edit

By now, we are sure that NDSolve has problems with DiracDelta. So there is another way to deal with it by approximating it by a NormalDistribution,

e1=0.00001;
p1 = NDSolve[{f''[x] == 
    PDF[NormalDistribution[10, e1], x] zf/(ep a^2) - 
     1 Sum[z[i] cp[i], {i, 1, 2}], f[10] == zf/(ep a^2), f'[20] == 0},
   f, {x, 10, 100}, MaxStepSize -> e1, 
  MaxSteps -> Infinity]

Finally, plotting the two results combine,

Show[Plot[f[x] /. p, {x, 10, 100}], 
 Plot[f[x] /. p1, {x, 10, 100}, PlotStyle -> {Dashed, Red}]]

enter image description here

Apparently, both numerical solutions p and p1 are identical.

$\endgroup$
  • $\begingroup$ @ MMM I am getting this error NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 10.?? $\endgroup$ – amin bk Dec 31 '16 at 19:46
  • $\begingroup$ @aminbk which part of the code is giving you this error? $\endgroup$ – zhk Jan 1 '17 at 3:45
  • $\begingroup$ @ MMM when I execute the code above I get the following error General::unfl: Underflow occurred in computation.General::stop: "Further output of \!(* StyleBox[ RowBox[{\"General\", \"::\", \"unfl\"}], \"MessageName\"]) will be suppressed during this calculation" $\endgroup$ – amin bk Jan 1 '17 at 12:15
  • $\begingroup$ when I execute p1 = .... do you know what is the reason? $\endgroup$ – amin bk Jan 1 '17 at 12:17
  • $\begingroup$ @aminbk Sorry dear. I have no clue. What version of MMA you are using? I have V11.0. $\endgroup$ – zhk Jan 1 '17 at 13:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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