# How to remove a singularity with WhenEvent

Bellow I have a differential equation which hits a singularity at low values of t. What I want to do is somehow utilize the WhenEvent command in order to replace the last factor of the equation (2/t a'[t]) with 0 at low values of t in order to avoid getting errors. Is there a way to do this?

Here is an example of the equation and how I attempted to use the WhenEvent Command.

q = NDSolve[{a''[t] ==
1/2 a[t]^2 + 1/6 a[t]^6 - 1/4 a[t]^4 - 2/t a'[t], a[0] == 1,
a'[0] == 0, WhenEvent[t < .1, a'[t] -> 0]}, a, {t, 0, 1}];


Any help would be appreciated!

-
Can you explain why you think you need a WhenEvent for this? Wouldn't e.g. this do what you need: NDSolve[{a''[t] == 1/2 a[t]^2 + 1/6 a[t]^6 - 1/4 a[t]^4 - Piecewise[{{2/t a'[t], t > 0.1}}], a[0] == 1, a'[0] == 0}, a, {t, 0, 1}] –  Albert Retey Apr 17 '13 at 16:27
I was thinking about an If[] instead. My guess for why WhenEvent is not working is because it works only at discrete events. Using an If[] or what @AlbertRetey should meet what you want –  Sosi Apr 17 '13 at 16:59
Please do not use the bugs tag when asking a new question. This tag is reserved for problems which have been confirmed by the community to be due to bugs in Mathematica. This does not appear to be the case here. –  Szabolcs Apr 18 '13 at 0:53
Sorry I think I accidentally tagged bug. Won't happen again. And thanks for your input! –  William John-Pierre Duhe Apr 22 '13 at 17:11

There are two problems with your code:

# 1. Incorrect Usage of WhenEvent

As the name indicates WhenEvent is meant to be used when you want to e.g. switch at certain events, what I think you try to do is to set a'[t] to zero for a whole period (0<=t<=0.1), but that's AFAIK not what WhenEvent can be used for directly. Of course you can reformulate your problem so that it will fit better for the usage of a WhenEvent, e.g. by introducing an artificial discrete variable which switches the singular term on only after t>=0.1:

NDSolve[{
a''[t] == 1/2 a[t]^2 + 1/6 a[t]^6 - 1/4 a[t]^4 - 2*k[t]/t* a'[t],
a[0] == 1, a'[0] == 0, k[0] == 0,
WhenEvent[t == 0.1, k[t] -> 1]
},
a, {t, 0, 1}, DiscreteVariables -> k
]


unfortunately this will not solve your problem: you'll still get Power::infy messages and nonnumeric intermediate results which make NDSolve fail. This is because of

# 2. 0/0 evaluates to nonumeric Indeterminate

which can be checked by evaluating 0/0. But as with the above formulation the whole term is switched off for 0<=t<0.1, we can now savely change it to something nonsingular for t<0.1, which will not change the result if the term gives the same results for t>=0.1. That will solve that problem, too:

NDSolve[{
a''[t] == 1/2 a[t]^2 + 1/6 a[t]^6 - 1/4 a[t]^4 - 2*k[t]/Max[t, 0.1]* a'[t],
a[0] == 1, a'[0] == 0, k[0] == 0,
WhenEvent[t == 0.1, k[t] -> 1]
},
a, {t, 0, 1}, DiscreteVariables -> k
]


You can check that this approach gives the same result as my suggestion from the comment (which of course could also be formulated with an If, When, Unitstep and a whole lot of other functions):

NDSolve[{
a''[t] == 1/2 a[t]^2 + 1/6 a[t]^6 - 1/4 a[t]^4 - Piecewise[{{2/t a'[t], t > 0.1}}],
a[0] == 1, a'[0] == 0
},
a, {t, 0, 1}
]


# Final Note

If you are interested, you can check that your original formulation will only trigger an event exactly once, at t=0.1. So it wouldn't help even if the 0/0 term would evaluate to zero with no errors. You can see what happens with the following code (probably you'll need to adjust the plot ranges...)

q4 = NDSolveValue[{a''[t] ==
1/2 a[t]^2 + 1/6 a[t]^6 - 1/4 a[t]^4 - 2/Max[t, 10^-10] a'[t],
a[0] == 1, a'[0] == 0}, a, {t, 0, 1}]

q3 = NDSolveValue[{a''[t] ==
1/2 a[t]^2 + 1/6 a[t]^6 - 1/4 a[t]^4 - 2/Max[t, 10^-10] a'[t],
a[0] == 1, a'[0] == 0, WhenEvent[t < .1,Print[t]; a'[t] -> 0]}, a, {t, 0, 1}]

Plot[{q3[t], q4[t]}, {t, 0, 1},
PlotRange -> {{0.08, 0.12}, {0.9999, 1.001}}]


-
Thanks so much for this help its fantastic! I have managed to tweak out some of the suggestions to get my program working! I will post a working model soon so you guys can see what I accomplished. –  William John-Pierre Duhe Apr 22 '13 at 17:09

I usually cope with such kind of singular problems in a different way, simple and practical, but not with constructs of the sort of the WhenEvent (which is anyway not for such cases, as it is shown in the Albert Retey's answer). I simply regularize the singularity. This, of course, is only possible in case it is not forbidden from the point of view of physics behind the problem being solved.

In the case of the equation above regularization means the replacement of 1/t by 1/Sqrt[t^2+eps], where eps is a small number, like 0.001. Its value may be chosen depending upon the proximity to zero your problem requires.

To illustrate it let us solve this system with the initial condition close to zero. Evaluate this:

q = NDSolve[{a''[t] ==
a[t]^2/2 + a[t]^6/6 - a[t]^4/4 - 2 a'[t]/(Sqrt[t^2 + 0.0001]),
a[0] == 0.01, a'[0] == 0}, a, {t, 0, 5}] // First;

Plot[Evaluate[a[t] /. q], {t, 0, 5},
AxesLabel -> {Style["t", Italic, 16], Style["a", Italic, 16]}]


It returns this:

-
Gee Alexei, don't leave us hanging... –  Ｊ. Ｍ. Apr 19 '13 at 8:32
@J.M. This your comment I do not understand. –  Alexei Boulbitch Apr 22 '13 at 7:10
"It returns this:" - what, precisely does the plot return in your machine? It certainly returns something different from the nice image you uploaded. –  Ｊ. Ｍ. Apr 22 '13 at 9:20
@J.M. It certainly returns what I uploaded, though I cannot see my uploaded images here. You may make it sure by evaluation the posted code. Is it your regular way to accuse unknown people in a lie? –  Alexei Boulbitch Apr 23 '13 at 12:18
I'm not saying you were lying; since you seem to very readily assume negative intent out of a simple nudge, let me speak out very directly as a moderator instead of as a regular user: I was prompting you to edit your post since we're not seeing what you intended us to see. –  Ｊ. Ｍ. Apr 23 '13 at 12:57