清空筛选条件
清空筛选条件

Why does this function appear noisy?

69 次查看（过去 30 天）
Ryan
Ryan 2024-9-6，10:42
编辑： John D'Errico 2024-9-6，16:38
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

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2024-9-6，12:22
编辑：John D'Errico 2024-9-6，12:28
Ran in:
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.
  2 个评论
Ryan
Ryan 2024-9-6，12:36
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
John D'Errico 2024-9-6，15:30
编辑：John D'Errico 2024-9-6，16:38
Ran in:
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Common Operations 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by