Why ''int'' and ''integral'' return different answers
10 ビュー (過去 30 日間)
Hello, I am doing a simple integration:
zzz = @(x) x.^(-0.98).*(1-x).^(-0.98)
While for int:
zzz = x.^(-0.98).*(1-x).^(-0.98)
I know 99.9361 is the right numer, but why integral returns a different number? (I prefer to use 'integral')
A big thank you in advance!
John D'Errico 2022 年 9 月 10 日
編集済み: John D'Errico 2022 年 9 月 10 日
As others have said, the problem with integral is because of the singularities. Numerical methods don't like singularities.
zzz = @(x) x.^(-0.98).*(1-x).^(-0.98);
So it is a good idea to avoid integral on functions like that. Unless, of course, you see that the integral you are trying to evaluate is just the complete beta function. It helps to know your special functions!
As such, you should see the integral is given very nicely by beta. No explicit integration was necessary. Don't forget to add 1 to the exponents to give as arguments to the beta function.
Alternatively, I could probably do the same computation using a Gauss-Jacobi quadrature, but that seems to be a bit of overkill for this problem, since beta already exists in MATLAB, and does what you want very nicely.
その他の回答 (3 件)
Torsten 2022 年 9 月 10 日
Your function has singularities at x = 0 and x = 1.
Integral doesn't manage to estimate the area under the curve (which goes to Infinity at both places) correctly.
Walter Roberson 2022 年 9 月 10 日
You should pass in abstol and reltol options to integral()
integral() subdivides intervals until the contribution appears to be small. But there are functions for which a small contribution adds up a lot.
Consider for example sum of 1/x as x approaches infinity, or the similar integral. For any given threshold you can find x such that all further 1/x will be less than the threshold. Let the threshold be eps(1) and you have the situation where if you were looping adding 1/n that each 1/n contribution would individually be less than a 1 bit difference in the floating point sum, so it would be reasonable to stop adding them. But integral 1/x is log(x) and as x approaches infinity that would go to log(infinity) which is infinite. This shows that numeric methods can be infinitely wrong if they do not happen to do the right transform to match the function.
It is not really a bug in the numeric integration, it is a limitation: for any numeric integration algorithm you can create functions that the algorithm cannot integrate accurately.
Bruno Luong 2022 年 9 月 10 日
編集済み: Bruno Luong 2022 年 9 月 10 日
I think the problem is more subttle than the simple presence of singularities. Some "mild" singularity can be handlle correctly by integration. The problem here is that for power that is close to -1, where the integration goes to infinity and it is not mild at all. It's more like of making a transition to diverging integration. No wonder why integrationn function cannot deal it well.
Normally for all w < 1, singularity exists but one observe that integration starts to have difficulty for w < 0.15, otherwise it converges and gives an accurate answer.
w = linspace(0.01,1);