Weird numerical integration behavior

1 回表示 (過去 30 日間)
Hussein Ammar
Hussein Ammar 2019 年 9 月 24 日
コメント済み: John D'Errico 2019 年 9 月 24 日
Hello all,
Please check this simple example:
d1 = 8;
d2 = 1.4142;
apPDF = @(x1) (1./sqrt(2*pi)) .* exp(-power((d2.*x1./sqrt(2)) - d1,2) ./ d2^2);
integral(apPDF, 0, Inf)
myVarBoundary = 10^24;
integral(apPDF, 0, myVarBoundary)
My Matlab output:
>> d1 = 8;
d2 = 1.4142;
apPDF = @(x1) (1./sqrt(2*pi)) .* exp(-power((d2.*x1./sqrt(2)) - d1,2) ./ d2^2);
integral(apPDF, 0, Inf)
myVarBoundary = 10^24;
integral(apPDF, 0, myVarBoundary)
ans =
1.0000
ans =
0
The correct answer that I should get from the second integration should be one, right?
What is wrong in my integration? Should I replace myVarBoundary with Inf after some value.
Best Regards,

採用された回答

John D'Errico
John D'Errico 2019 年 9 月 24 日
編集済み: John D'Errico 2019 年 9 月 24 日
No. You should NOT make it inf. Even 1e24 is incredibly large, wild overkill.
As to why you are using integral to compute the area under what loos like a normal PDF is completely beyond me.
d1 = 8;
d2 = 1.4142;
apPDF = @(x1) (1./sqrt(2*pi)) .* exp(-power((d2.*x1./sqrt(2)) - d1,2) ./ d2^2);
You need to understand what integral does when it sees a function with limits that wide. It evaluates the function at a variety of points in the interval. Lets try a few, just for kicks.
apPDF(0)
ans =
5.04917109360107e-15
>> apPDF(1e24)
ans =
0
>> apPDF(1e24 / 2)
ans =
0
>> apPDF(1e24 / 100000)
ans =
0
>> apPDF(1e24 / 10000000000000)
ans =
0
Do you see anything significant? Do you see a function that seems to be everywhere zero on the interval [0,1e24]? And even when not identically zero, it deviats from zero on the order of the convergence tolerance. Should it somehow, magically know that in effectively a tiny corner of that HUGE interval, it is non-zero?
apPDF(5)
ans =
0.00443082846737379
So now if we do this:
integral(apPDF,0,100)
ans =
1
What a surprise! It integrates to 1.
  3 件のコメント
Hussein Ammar
Hussein Ammar 2019 年 9 月 24 日
Hello all, thank you for your replies.
I understand what your pointing for, but unfortunately I can't control the value of "myVarBoundary" (depending on a lot of things it can be very small or very large). so, I think I can assume that if "myVarBoundary" crosses some threshold, my answer should tend to 1. right?
Anyway, I think what's wrong is to use myVarBoundary directly as the upper limit of the numerical integration, because when myVarBoundary is very large, the numerical integration won't see the very tiny part that is non-zero at the start of integration domain.
Regards,
John D'Errico
John D'Errico 2019 年 9 月 24 日
Steve makes an excellent point. Here, you might need to be looking at the ground using the Hubble space telscope though, to get the necessary resolution.
As I showed, the function was zero above x1=100. So out of an interval of width 1e24, it is zero on only the fraction (well) below 100. 1 part in 1e22?
1e22 is a number almost as large as Avogadro's number. Big.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeDebugging and Analysis についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by