Why does this function appear noisy?

Hello,

I'm troubleshooting a numerical integration and have narrowed the problem down to a particular part of my integral kernel, consider the following functions f1 and f2 and parameters for setup

p = 10^12;
gm1 = 100;
tfinal = (10^2)/p;
tsize = 50;
c= 3*10^8;

beta1 = sqrt(1-1/gm1^2);
v = 2.5*10^3; 
t = tfinal/tsize:tfinal/tsize:tfinal;
gm = sqrt(1+(c/v*(t-t(1))+gm1*beta1).^2);
beta = sqrt(1-1./gm.^2);
f1 = @(s,sp,k) ((s-sp)/v+gm(k)-beta(k).*gm(k)).*((s-sp)/v+gm(k)+beta(k).*gm(k));
f2 = @(r,s,sp,k) (1+r.^2/v^2+f1(s,sp,k));

testFunc = @(r,s,sp,k) f2(r,s,sp,k).^2 -4*f1(s,sp,k)

Using fplot produces the following noisy result, which worsens with increasing "v"

fig1=figure; fplot(@(sp) testFunc(0,0,sp,10), [-5*10^-6 5*10^-6])

Conversely, we can show analytically that for r=0 testFunc is equivalent to

testFunc2 = @(s,sp,k) (f1(s,sp,k)-1).^2
fig1=figure; fplot(@(sp) testFunc2(0,sp,10), [-5*10^-6 5*10^-6])

which I believe is evidently what the function should look like for r=0 plotted along sp.

Presumably the source of the noise is some issue with numerical precision, but I do not know precisely what is the cause or how to address it.

Thanks

Answer 1

[John D'Errico]

Do you understand how floating point arithmetic works? That you only have a finite number of digits, and that when you subtract two numbers that are almost equal, then you get in trouble? MATLAB does not compute in infinite precision. This is commonly called massive subtractive cancellation. For example, if we compute the difference of two numbers, perhaps

a1 = @(x) 1 + x;
a2 = @(x) 1 + sin(x);

Now, when x is very small, we might expect to see a simple, smooth curve. After all, both functions are well defined, right?

fplot(@(x) a1(x) - a2(x),[-1e-5,1e-5])

The difference between the two functions should look almost live a cubic polynomial. What happened? Massive subtractive cancellation. Note that to make things really bad, I added both x and sin(x) to 1. Instead, plot this on the same interval.

fplot(@(x) x - sin(x),[-1e-5,1e-5])

Mathematically, the two expressions to plot are identical. But in terms of floating point arithmetic, they are not even close. I would note that on a narrower interval yet, even forming the differece of x and sin(x) will cause problems.

fplot(@(x) x - sin(x),[-1e-7,1e-7])

What did you do at the very end?

% testFunc = @(r,s,sp,k) f2(r,s,sp,k).^2 -4*f1(s,sp,k);

So you created the DIFFERENCE of two expressions that are very nearly equal to each other. The result is expected to be quite small.

% fig1=figure; fplot(@(sp) testFunc(0,0,sp,10), [-5*10^-6 5*10^-6])

And you saw floating point trash in the result. SURPRISE!

What could you have done? You might consider using symbolic tools to perform the same computations. Now you could work in sufficiently high precision to avoid the floating point trash.

[Ryan]

Hello, thanks for your response, I see my error and have since found that rewriting the expression eases things slightly.

Does the symbolic toolbox work to a higher precison then?

[John D'Errico]

Yes. Well, you can make it do so, using tools like digits and vpa, where you can control the precision. Don't forget, it will also be WAY slower. Or you could use my HPF toolbox, which would also be really relatively slwo compared to doubles. The price you pay for higher precision is in computing time. You also need to learn to be careful about mixing doubles and syms, as sometimes you get strange stuff.

For example, the problem I showed before is now trivially done using syms. I never even needed to force it to do anything with explicitly higher precision.

syms x
a1(x) =  1 + x;
a2(x) =  1 + sin(x);
fplot(a1(x) - a2(x),[-1e-5,1e-5])

How about your problem?

p = sym(10^12);
gm1 = sym(100);
tfinal = (10^2)/p;
c = sym(3*10^8);
beta1 = sqrt(1-1/gm1^2);
v = sym(2.5*10^3); 

All of that is trivial.

syms t
syms s sp r
tstart = -tfinal;
gm = sqrt(1+(c/v*(t-tstart)+gm1*beta1).^2);
beta = sqrt(1-1./gm.^2);

f1 = ((s-sp)/v+gm-beta.*gm).*((s-sp)/v+gm+beta.*gm);
f2 = 1+r.^2/v^2+f1

f2 = 

testFunc = f2.^2 -4*f1

testFunc = 

So far, I did nothing beyond pure symbolics. I'll plug in the values of r and s next.

TF_sp = subs(testFunc,[r,s],[0,0])

TF_sp = 

Note this is still a function if sp, and of t. In your case, when k==10, then t(10) was 2e-11. So far, there has been really nothing done in terms of floating point arithmetic. YET. Next, we can turn this mess into something with floating point numbers in it. I'll just show it using only 5 digits on the coefficients so we can see it for what it is.

vpa(expand(subs(TF_sp,t,sym(2e-11))),5)

ans = 

In fact, the result for a fixed value of t was simply a 4th degree polynomial. The problem is, fplot decides to use doubles to evaluate it.

fplot(subs(TF_sp,t,sym(2e-11)),[-5*10^-6 5*10^-6])

No matter what it does though, that polynomial has coefficients that vary HUGELY. And that will be a numerical problem.

Can we do better? Perhaps. One trick might be to use a transformation of variables. I'll use a substitution and see what happens to the coefficients.

 sphat = sp/1e-6;

Now we would have the polynomial

vpa(expand(subs(TF_sp,[sp,t],[sphat*1e6 , sym(2e-11)])),5)

ans = 

Do you see now what is happening? When sp is on the order of 1e-6, then sphat will be on the order of 1. We will need to carry dozens of digits to evaluate that mess. At least 50 digits, and probably more to be safe.

Slightly better is to use a different transformation, one that is intermediate, chosen to make the coefficients more tractable.

sphat = sp/1e-3;
vpa(expand(subs(TF_sp,[sp,t],[sphat*1e3 , sym(2e-11)])),5)

ans = 

Do you see I chose one to make the higher order coefficients all have roughly similar magnitudes.

fplot(subs(TF_sp,[sp,t],[sphat*1e3,sym(2e-11)]),[-.005,.005])

Now the roughness goes away, and I never needed to use super high precision.

Ok, I could probably have gotten it to work in a different way too.