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?

Improve this question
$\endgroup$
3
  • 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
    Commented Dec 31, 2016 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
    Commented Dec 31, 2016 at 14:40
  • $\begingroup$ Possible duplicate: (55986). Also related: (133493) $\endgroup$
    – Michael E2
    Commented Jan 1, 2017 at 4:20

3 Answers 3

Reset to default
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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks alot, yes , the problem for me is why mathematica does not understand this. $\endgroup$
    – amin bk
    Commented Dec 31, 2016 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
    Commented Dec 31, 2016 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

Improve this answer
$\endgroup$
2
  • $\begingroup$ Why the solution curve starts in the second quadrant? $\endgroup$
    – zhk
    Commented Dec 31, 2016 at 16:07
  • 1
    $\begingroup$ @MMM. Problem is a Plot gives strange output.I'm edited answer. $\endgroup$
    – Mariusz Iwaniuk
    Commented Dec 31, 2016 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.

Improve this answer
$\endgroup$
8
  • $\begingroup$ @ MMM I am getting this error NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 10.?? $\endgroup$
    – amin bk
    Commented Dec 31, 2016 at 19:46
  • $\begingroup$ @aminbk which part of the code is giving you this error? $\endgroup$
    – zhk
    Commented Jan 1, 2017 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
    Commented Jan 1, 2017 at 12:15
  • $\begingroup$ when I execute p1 = .... do you know what is the reason? $\endgroup$
    – amin bk
    Commented Jan 1, 2017 at 12:17
  • $\begingroup$ @aminbk Sorry dear. I have no clue. What version of MMA you are using? I have V11.0. $\endgroup$
    – zhk
    Commented Jan 1, 2017 at 13:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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