My Math Forum  

Go Back   My Math Forum > Math Forums > Math Software

Math Software Math Software - Mathematica, Matlab, Calculators, Graphing Software

LinkBack Thread Tools Display Modes
March 22nd, 2017, 02:44 PM   #11
Senior Member
Joined: Jul 2008
From: Western Canada

Posts: 2,951
Thanks: 33

I'll have a look at Pari. In fact, I think I already have it on my computer. I recall downloading it quite a while back.

Originally Posted by studiot View Post
Have you tried the Second Euler Summation Formula (used for improper integrals, particularly trigonometric ones) ?...Press and Flannery
I have a copy of Numerical Recipes in C. In the section on improper integrals, they give the Second Euler-Maclaurin summation formula which I assume is the one you're referring to. They don't mention that it's good for trigonometric functions. So I overlooked it. Also, I utterly despise the way that they've written the code, with an apparent attempt to obscure what they've done, so that you have to buy a license to use their version rather attempt to rewrite your own. However, I'll set aside my bias, and see if I can get it to work for me.<\end rant>
Thanks for pointing it out.
Yooklid is offline  
March 26th, 2017, 11:47 PM   #12
Senior Member
Joined: Jul 2008
From: Western Canada

Posts: 2,951
Thanks: 33

billymac00, There seems to be something corrupted in my user settings, which makes private messages unpredictable. So, I thought it best to mention that I did get your PM, and hopefully you got my reply. If not, just let me know here in this thread.

I haven't yet tried your Pari code, because I'll have to figure out where I actually installed Pari, and also because I decided to pursue studiot's suggestion of the Second Euler Summation Formula. This led to other things, which in turn, led to even more things (branches going off in all directions), all of which turned out to be very interesting.

Regarding Euler-Maclaurin, I did some research. Looking at an article on the Euler-Maclaurin formula, I had one of those epiphany moments, and realized that—for a completely unrelated problem that I'm also working on—this would give the solution for a finite summation that I'd been trying to solve. So, I solved one problem, but not the one that I originally posed here.

On the subject of the Second Euler Summation Formula for the original problem, from what I had read, there didn't appear to be a huge difference between this and the basic trapezoidal rule. Further reading of some interesting articles on the trapezoidal rule indicate that for cyclical functions, the trapezoidal rule may be optimum, giving exponential convergence. I thought I'd tried the trapezoidal rule back in the early days when I first worked on this problem. But I decided to revisit it. I confirmed that it didn't converge nearly as fast as my customized Simpson's rule solution, and I wondered why.

The reason why the trapezoidal rule is supposed to converge so fast for cyclical functions is that the approximation error cancels out more quickly than the worst case error formulae would predict. The trapezoidal rule underestimates where the function is concave downwards, and overestimates where the function is concave upwards. Since a cyclical function alternates between concave up and down, much of the error cancels, giving very fast convergence. But curiously, not in my case.

I decided that I should have another look at the function I'm trying to integrate. And on that note, I'll point out that I made a typo when I posted the function. I had a $cos^2$ term in the numerator, and it should have been just a first order $cos$. This is the corrected expression:

$f(x)=\dfrac{k_1 cos(x-k_2)+k_3}{\sqrt{k_4(k_5-k_5 cos(x-k_2))+k_6x}}$

When I first started working on this problem about 5 years ago, I graphed the function for a particular set of $k_n$ values, and nothing appeared amiss. It looked like a typical decaying cosine function. However, I've now realized that with certain values of $k_n$, this can change so that the graph looks like a series of extremely narrow spikes occurring at intervals of $2π$. This explained why the concave up and concave down won't cancel. I realized that to solve this efficiently, I need to concentrate the function samples around multiples of $2π$.

I decided to split the integral into sub-integrals with limits $0..\pi, \pi..3\pi, 3\pi..5\pi, 5\pi..7\pi$, etc., so that the spike always occurs in the centre of the integration limits (except for the first, and possibly last interval). Then I applied the double exponential method for each sub-integral (which from my past experience has been most efficient for these kinds of integrals). The result is that I'm seeing a reduction in the required number of function evaluations by a factor of about 10 or more to achieve the same error. Considering what the graph of the function looks like, I suspect I won't get much better than this.

Also of note (since I'm not integrating off to infinity, fortunately): As can be seen from the formula above, the spikes decay very slowly as a function of $1/\sqrt{x}$, which means that I have to count every one, which is 400 in the case that I'm currently working on.

Anyway, I'd like to say that I appreciate the suggestions you've all given me. They may not necessarily have led to what you expected, but it's had very positive results. I still intend to pursue some of the suggestions that I haven't yet followed. I'm already thinking of a couple tangents to veer off on.
Yooklid is offline  

  My Math Forum > Math Forums > Math Software

algorithm, integration, numerical

Thread Tools
Display Modes

Similar Threads
Thread Thread Starter Forum Replies Last Post
Help with Numerical Integration Please :) AnomanderRakeSoD Applied Math 0 June 21st, 2013 02:11 AM
numerical integration Albert.Teng Calculus 3 July 1st, 2012 10:13 PM
Numerical integration in Excel dallairius Calculus 1 April 6th, 2012 09:07 AM
Question on numerical integration PageUp Applied Math 1 March 19th, 2012 05:58 AM
[Numerical Integration] Change of variable VFernandes Applied Math 1 January 5th, 2012 12:38 PM

Copyright © 2017 My Math Forum. All rights reserved.