1
$\begingroup$

I have the following function

 f[s_]:=(11/10 + s) (-((50 (-(13/10) + (169 s)/100))/( 13 Sqrt[231] Sqrt[ 100/33 + s] (1/4 + (11 s)/10 + s^2))) - 1/( 4 (11/10 + s)^( 3/2) (s + 1/(4 (11/10 + s)))))^2

Which is continuous and well-defined for

-1<=s<=0

However, at s=sp, there is a very non-trivial cancellation of singularities, where :

sp=1/2 (-(11/10) + Sqrt[21]/10)

This is a problem when manipulating the function numerically, as can be seen by plotting it for sp-0.0001<s<sp+0.0001, where we see huge oscillation due (I assume) to the rounding errors of cancelling two big factors.

This wouldn't be a problem, but I am using this function as a distribution density, and when using RandomVariate to generate points, it will slow down immensely when drawing close to sp. For this reason, I would like to smooth out this function around sp. Would there be an elegant way to do it ?

The only thing I could come up with is define a piecewise function and interpolate linearly between sp-0.0001 and sp+0.0001, but I was wondering if there was a better way to do this, maybe with a built-in function that I ignore.

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ I tried, but that doesn't help. Indeed, the cancelling of the divergence is not a simple one such as $x^2/x$. While the two pieces of the functions diverge in the same and opposite way at $sp$, these factors differ away from sp. In other words, only the first (singular) term in an expansion around $sp$ cancels, but this does not hold for all s. (Original message was deleted, it was suggesting the use of "Together". I will leave my comment as I think it is an important clarification). $\endgroup$ Commented Oct 17, 2022 at 15:29
  • $\begingroup$ Your expression is equal to (1690 (20-13 s)^2)/(100+33 s)*1/(20 Sqrt[77+70 s]-26 s Sqrt[77+70 s]+7 Sqrt[1000+330 s])^2 if I did not make a mistake. In this way of writing it, there is no cancellation between -1 and 0. $\endgroup$ Commented Oct 17, 2022 at 15:35
  • $\begingroup$ Indeed, that seems to be the case. May I ask how you arrived at this value ? Indeed my function is actually dependent on several parameters, so I would need to do that systematically (Edit : I found it, it is simply simplify. I will try that out if it works generally) $\endgroup$ Commented Oct 17, 2022 at 15:40

2 Answers 2

Reset to default
3
$\begingroup$ 
Clear["Global`*"] f[s_] := (11/10 + s) (-((50 (-(13/10) + (169 s)/100))/(13 Sqrt[231] Sqrt[ 100/33 + s] (1/4 + (11 s)/10 + s^2))) - 1/(4 (11/10 + s)^(3/2) (s + 1/(4 (11/10 + s)))))^2 // Simplify sp = 1/2 (-(11/10) + Sqrt[21]/10); f[sp] (* ComplexInfinity *)

Define the value at sp as the Limit

f[sp] = Limit[f[s], s -> sp] // FullSimplify (* -((507 (-27949795 + 2346512 Sqrt[21]))/138358705156) *)

To Plot in the vicinity of sp, specify a WorkingPrecision to avoid machine precision calculations and increase MaxRecursion

Plot[f[s], {s, sp - 0.0001, sp + 0.0001}, WorkingPrecision -> 15, MaxRecursion -> 5]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you for the trick of defining at $sp$ with the Limit. However, the workaround with the workingprecision and recursion is for the Plot. However, I would like to smooth out the function for the purposes of applying RandomVariate on it efficiently. I should have probably been more precise in the question, I will try to amend it if possible. ( I will try just defining well on sp, wether that improves the randomVariate) $\endgroup$ Commented Oct 17, 2022 at 15:43
  • $\begingroup$ You can also specify a WorkingPrecision with RandomVariate and RandomReal $\endgroup$ Commented Oct 17, 2022 at 15:45
  • $\begingroup$ Will try that ! Edit : Adding explicitely the WorkingPrecision in RandomVariate (although I though that 15 was the default) indeed solves also the problem and evaluates fast. I will mark this question as accepted as it applies more broadly to the problem, unlike the other valid answer which is specific to my function. $\endgroup$ Commented Oct 17, 2022 at 15:58
2
$\begingroup$

Start with OPs expression

f1 = (11/10+s) (-((50 (-(13/10)+(169 s)/100))/(13 Sqrt[231] Sqrt[100/33+s] (1/4+(11 s)/10+s^2)))-1/(4 (11/10+s)^(3/2) (s+1/(4 (11/10+s)))))^2;

Simplify

f2=f1//FullSimplify

enter image description here

The problem is that at the point that OP mentions, both the numerator factor and the factor $5+22s+20s^2$ vanish. Define the auxiliary

A=20 Sqrt[77+70 s]-26 s Sqrt[77+70 s]+7 Sqrt[1000+330 s];

This is like the numerator factor except for some signs that I changed (I tried different possibilities). The idea is to write f2 as (f2*A^2)*(1/A^2) and to hope that f2*A^2 does not have those singularities. In fact

g=f2*A^2//FullSimplify

enter image description here

Therefore f1 is equal to

g*(1/A^2)

enter image description here

but in this way of writing it, the singularities are absent.

Improve this answer
$\endgroup$
3
  • $\begingroup$ This seems to work, it appears that the singularities can actually cancel out, contrary to my belief. For some reason, simply using Simplify and using the function as indicated in the first step also works, although I don't know why. $\endgroup$ Commented Oct 17, 2022 at 15:53
  • $\begingroup$ Are you saying that with f1//Simplify you also get the expression that I wrote at the end? (I do not.) If not, what do you mean? $\endgroup$ Commented Oct 17, 2022 at 15:55
  • 1
    $\begingroup$ I do not (I get the one you wrote at the start), but plotting the function close to $sp$ yields no numerical artifacts. $\endgroup$ Commented Oct 17, 2022 at 16:00

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.