Threaded Mode | Linear Mode

ppphzj · 08-27-2024, 12:43 PM

the 3924 is wrong

the 3844 is right

Nigel (UK) · 09-02-2024, 06:44 PM

The different results given by v3844 and v3924 are due to the change in the integration algorithm used, which happened in v3893. Before that the Romberg method was used; the current code uses the Double Exponential method. (If you wish you can read about this here.)

For many integrals the new method is quicker and more accurate. It is particularly good when there is an integrable singularity at either end of the interval of integration. However it is less good when integrating functions that vary periodically many times over the interval, which is the case for your integral.

Of course, since your integrand is periodic it’s easy to evaluate it by integrating over one cycle, multiplying by the number of complete cycles, and then adding the result of integrating over the final partial cycle. In general this won’t be possible and you will have to be content with spotting the sorts of integrand that are likely to give problems.

I don’t know an integration method which is superior to all other methods for all types of integrand. The double exponential method is good but - as your example shows - not perfect. It would be possible to re-include the Romberg code in addition to the new code - what do you think? How much code space would it be worth sacrificing for this?

Nigel (UK)

ppphzj · 09-03-2024, 10:45 AM

(09-02-2024 06:44 PM)Nigel (UK) Wrote: The different results given by v3844 and v3924 are due to the change in the integration algorithm used, which happened in v3893. Before that the Romberg method was used; the current code uses the Double Exponential method. (If you wish you can read about this here.)

For many integrals the new method is quicker and more accurate. It is particularly good when there is an integrable singularity at either end of the interval of integration. However it is less good when integrating functions that vary periodically many times over the interval, which is the case for your integral.

Of course, since your integrand is periodic it’s easy to evaluate it by integrating over one cycle, multiplying by the number of complete cycles, and then adding the result of integrating over the final partial cycle. In general this won’t be possible and you will have to be content with spotting the sorts of integrand that are likely to give problems.

I don’t know an integration method which is superior to all other methods for all types of integrand. The double exponential method is good but - as your example shows - not perfect. It would be possible to re-include the Romberg code in addition to the new code - what do you think? How much code space would it be worth sacrificing for this?

Nigel (UK)

Thank you very much for your explanation. The new version is really much faster than the 3844 version, I will use the faster one.

Albert Chan · 11-13-2024, 12:19 AM

(09-02-2024 06:44 PM)Nigel (UK) Wrote: For many integrals the new method is quicker and more accurate. It is particularly good when there is an integrable singularity at either end of the interval of integration. However it is less good when integrating functions that vary periodically many times over the interval, which is the case for your integral.

You are right. DE algorithm does not like oscillatory integrand.
It assumed transformed function is bell-shaped, then simply sum the rectangles.

For this case, DE algorithm work too. It is just that software had placed a hard limit on levels deep.
Note: my code changed meaning of levels, start from 1 instead of 0 --> levels = 12 get 12 rows output.

lua> f = fn'x: exp(sin(x)*cos(x))'
lua> ; Q.quad(f, 0, 1000, 1, 12, 1e-16, true)

Code:

4366650221054      0.370273314848271

093244328251       0.055352241458539

2599753713897       0.17247834220460842

9010969740907      0.027945254119896016

6160848617153      0.0071589003168768695

4116213556815      0.0394516500021799

105792170186       0.004080945453478147

9227093485367      0.01557679187876928

9411055724377      0.02259616840292012

8404586620848      0.0008453834241068643

8468942997813      6.049402156368027e-06

8468944036285      9.76149629141775e-11

Problem is this take about 2^(levels + 2) = 2^14 samples.

Still, it beat Romberg by a factor of 2

lua> ; Q.integ(Q.simpson(Q.u(f, 0, 1000)), 2^15)

Code:

4       865.1840001860048       299.21141073191774

8       778.0020184493249       693.544473641916

16      896.6700594243979       1014.4110063293531

32      1107.3997862356239      1317.7179315491724

64      1035.9832383865798      964.6015618987897

128     1049.6893911940772      1063.3938708872174

256     1084.9625972289223      1120.234726810947

512     1057.8616940259467      1030.7609975864543

1024    1064.5204424618387      1071.1791781971756

2048    1063.8045493118993      1063.0886565033247

4096    1063.8480952794305      1063.8916412417711

8192    1063.8468965924299      1063.8456979054642

16384   1063.846894123814       1063.846891655198

32768   1063.8468944049175      1063.8468946860226

col 1 = subintervals
col 2 = extrapolated trapezoid areas
col 3 = extrapolated rectangle areas, height = midpoint
col 3 are less good, but it cost nothing to generate, and useful to cross check col 2