I want to solve the following expression (used to obtain an analytic solution to a current distribution inside a workpiece):

$$a_{mn} = -\frac{\frac{4}{ab} \int_0^a \int_0^b f(x',y')\sin(px')\sin(qy')\mathrm{d}x'\mathrm{d}y'}{t\sinh(tc)}$$

Here. $a$,$b$ are scalar constants and $p = \frac{m\pi}{a}$, $q = \frac{n\pi}{b}$, $t=\sqrt{p^2+q^2}$.

The function $f$ is a Gaussian distribution: $f(x,y) = \frac{I_0d}{\pi\sigma^2} \exp(-\frac{r^2d}{\sigma^2})$

I was wondering if using the Simpson's method is the smartest way to solve the double integral. Are there alternative solutions which would be more computationally efficient?

Additional information: The resulting matrix $a_{mn}$ is used in the following infinite series to obtain the final result: $\sum_{m=1}^\infty \sum_{n=1}^\infty a_{mn}p\cos(px)\sin(qy)\cosh(tz)$


It depends on how big $m$ and $n$ are (assuming that $x'$ ranges from $0$ to $a$, though the notation suggests the opposite if read literally). You should probably split the intervals into at least $m$ resp $n$ subintervals, integrate the pieces separately, and sum the results.

If $f$ is smooth, a Gaussian rule is probably more accurate for the same effort.

If $m$ and $n$ are large, it is better to rewrite the integral in terms of Fourier transforms and solve it by FFT.

| cite | improve this answer | |
  • $\begingroup$ Thanks! I had accidentally omitted the definition of f. I updated the question now. Indeed it is a smooth function, and indeed x does run from 0 to a. The result of this integral is fed into an infinite series, so in theoretically m,n are infinite (in practice I hoped to get by with a big number cut-off). $\endgroup$ – akid Jul 13 '12 at 12:24
  • $\begingroup$ @akid: see my edit. $\endgroup$ – Arnold Neumaier Jul 13 '12 at 14:36
  • $\begingroup$ Thanks again. To be honest, using a Fourier transform was an answer I was expecting. Unfortunately, I don't know much about them. Can you extend your answer with a note or a link on how to get started? $\endgroup$ – akid Jul 15 '12 at 12:46
  • $\begingroup$ en.wikipedia.org/wiki/Fourier_transform is a meaningful start. You really need to understand Fourier transforms if you work on the kind of problems you describe. The time spent on learning it is very well invested. It helps if you transform your trigonometric functions into exponentials and work in the complex domain - everything simplifies conceptually. $\endgroup$ – Arnold Neumaier Jul 15 '12 at 12:53
  • 4
    $\begingroup$ @akid: As you'll notice when you review your past experiences, whatever you do is at first a one-time thing, but spending the time to master it often makes it the beginning of a many-time thing. It is hardly ever in vain to learn properly even what one needs to do only once. - Knowing how to write trig functions as exponentials may even help you to rewrite your infiniteseries into a simpler form (interchange sum and integral). $\endgroup$ – Arnold Neumaier Jul 16 '12 at 7:23

If you solve this only for a single pair $m,n$, then high order Gaussian integration is likely the fastest way to do this. But if you want to build the elements of a matrix $a_{mn}$ for many different values $m,n$, then it is best to recognize that what you're computing here is a double Fourier transform, and it would likely be best to compute the double integral exploiting this fact using, for example, the FFT.

| cite | improve this answer | |
  • $\begingroup$ Wolfgang, a similar comment for you as for @Arnold: I did suspect a Fourier transformation, but I'm not not accustomed with them. Can you extend your answer with a link or hint how to get started? $\endgroup$ – akid Jul 15 '12 at 12:47
  • $\begingroup$ I think Arnold already answered this well. $\endgroup$ – Wolfgang Bangerth Jul 15 '12 at 15:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.