### Author:  [@Pier](https://discourse.julialang.org/u/Pier)
#### Post date: [May 16, 2019, 3:54pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/2 "2019-05-16T15:54:12Z")
Do you have linear/surface/volume current distributions? I have very recently done some Biot-Savart integrations for arcs/bars with different quadrature libraries in Julia and I may be able to provide assistance.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 16, 2019, 5:15pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/3 "2019-05-16T17:15:58Z")
> [@wiktorkujawa](#):
>
> Before using Julia i did it in C by using Complex Newton-Cotes formulas(Boole rule, Simpson rule and others):
Unless you are limited to equally spaced points for some reason, I would normally recommend using methods like Gaussian quadrature or Clenshaw–Curtis quadrature that use unequally spaced points and obtain much faster convergence for smooth integrands. Packages for this in Julia include
- [GitHub - JuliaMath/QuadGK.jl: adaptive 1d numerical Gauss–Kronrod integration in Julia](https://github.com/JuliaMath/QuadGK.jl) (adaptive)
- [GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for Gaussian quadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)
- [GitHub - JuliaMath/HCubature.jl: pure-Julia multidimensional h-adaptive integration](https://github.com/stevengj/HCubature.jl) (multidimensional, adaptive)
- [GitHub - giordano/Cuba.jl: Library for multidimensional numerical integration with four independent algorithms: Vegas, Suave, Divonne, and Cuhre.](https://github.com/giordano/Cuba.jl) (multidimensional, adaptive)
- [GitHub - robertdj/SparseGrids.jl: Sparse grid quadrature in Julia](https://github.com/robertdj/SparseGrids.jl) (multidimensional sparse grids)
- [GitHub - JuliaMath/Cubature.jl: One- and multi-dimensional adaptive integration routines for the Julia language](https://github.com/stevengj/Cubature.jl) (multidimensional and Clenshaw–Curtis rules, adaptive)
> [@wiktorkujawa](#):
>
> However I want to use GPU to speed up my programs.
Before you think about GPUs or parallelization, make sure to think about the underlying algorithms! Gaussian quadrature etcetera can be _exponentially_ faster than Simpson’s rule, for example.
Also, if you are evaluating Biot–Savart integrals over and over for the same current loop but different points in space, there are potentially _many_ additional ways to accelerate things, e.g. [fast multipole algorithms](https://ieeexplore.ieee.org/document/5754906) (also [here](https://ieeexplore.ieee.org/document/4202642)), interpolation schemes, semi-analytical methods to remove the singularity for points near the current source, etcetera.
---
### Author:  [@wiktorkujawa](https://discourse.julialang.org/u/wiktorkujawa)
#### Post date: [May 16, 2019, 7:10pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/4 "2019-05-16T19:10:03Z")
I deal with various current distributions. That is more technical(but in 3-dimension) like overhead lines, underground cables, utility poles and transmission towers on actual energetic infrastructure objects. So my calculations need to tell apart on Cartesian(X,Y,Z) coordinates. I even write application to calculate this, but that was console application and used only CPU in these calculations, but it works fine however i am sure that this can work much better.
---
### Author:  [@wiktorkujawa](https://discourse.julialang.org/u/wiktorkujawa)
#### Post date: [May 16, 2019, 7:33pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/5 "2019-05-16T19:33:36Z")
Thanks for help. If it goes about my problem I deal with various, nonsymmetrical problems and I wrote console and fine working program based on [this article](http://pe.org.pl/articles/2015/7/31.pdf) with using of this Boole Rule(in use in electro-energetic infrastructure) and even make a engineer thesis of whole question(calculation, mathematics and physics issues, application, environment influence of EM Field, engineering use). However I am sure that this can work better.
---
### Author:  [@Pier](https://discourse.julialang.org/u/Pier)
#### Post date: [May 17, 2019, 2:21pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/6 "2019-05-17T14:21:45Z")
I see from the paper that you compute the fields from wire configurations which require 1-D integrations. In this case I suggest to use the QuadQK package ([https://github.com/JuliaMath/QuadGK.jl](https://github.com/JuliaMath/QuadGK.jl)) which impements state-of-the art adaptive integration routines, with the benefit that it is very well suited to integrands with singularities (which can be passed to the routine). To increase performance I would not resort to GPUs but rather parallelize the computation within standard Julia by either parallelizing on the conductor pieces (for all points) or on the evaluation points (for all conductors). Hope this helps.
---
### Author:  [@wiktorkujawa](https://discourse.julialang.org/u/wiktorkujawa)
#### Post date: [May 17, 2019, 5:19pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/7 "2019-05-17T17:19:42Z")
I thinked to separate blocks for evaluation points and threads of the blocks for all segments of conductors.
---
### Author:  [@wiktorkujawa](https://discourse.julialang.org/u/wiktorkujawa)
#### Post date: [May 21, 2019, 8:30am UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/8 "2019-05-21T08:30:41Z")
I tested this package and time needed to compute integrals. I don’t know why, but it always(except first run) last about 0.9 s. However when I reduce rtol parameter to very small values it compile longer. My C algorithm last about tens milliseconds. They both give good results(too good so I can’t compare their precision). Maybe you know some integral that is hard too calculate on which I could test these algorithms.
The second question is about output of quadgk function, because I get two values result and error. From my point of view is good for testing, but if I want to use it in production I need only result value, so it’s waste of time, because algorithm have to do some unnecessary operations to get the second value - error.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [May 21, 2019, 11:35am UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/9 "2019-05-21T11:35:54Z")
> [@wiktorkujawa](#):
>
> I don’t know why, but it always(except first run) last about 0.9 s. However when I reduce rtol parameter to very small values it compile longer.
It is hard to say anything about this without a self-contained example.
> [@Please read: make it easier to help you](https://discourse.julialang.org/t/psa-make-it-easier-to-help-you/14757):
>
> Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …
> [@wiktorkujawa](#):
>
> I want to use it in production I need only result value, so it’s waste of time
It is a very cheap operation (using the previous refinement of the algorithm, if I remember correctly).
That said, it is prudent to use it in production too, eg error if the error is above some threshold that your calculations would tolerate.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 21, 2019, 11:55am UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/10 "2019-05-21T11:55:39Z")
> [@wiktorkujawa](#):
>
> The second question is about output of quadgk function, because I get two values result and error. From my point of view is good for testing, but if I want to use it in production I need only result value, so it’s waste of time, because algorithm have to do some unnecessary operations to get the second value - error.
The error estimate is produced automatically as a byproduct of the algorithm, with no additional function evaluations.
---
### Author:  [@wiktorkujawa](https://discourse.julialang.org/u/wiktorkujawa)
#### Post date: [May 21, 2019, 6:10pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/11 "2019-05-21T18:10:19Z")
Ok, I checked some example and now I know that is not that bad as I thinked.
I use vscode and when I compile this:
`@time integral, err = quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)`
I get:
```julia
@time integral, err = quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)
0.839247 seconds (1.42 M allocations: 86.103 MiB, 3.32% gc time)
(0.746824132812427, 7.887024366937112e-13)
```
but when I compile this:
```julia
@time for i=1:100
integral, err = quadgk(x -> exp(-x^2), 0, i, rtol=1e-8)
println("Result: $integral")
end
```
I get(I’ll shorten my output 🙂):
```julia
@time for i=1:100
integral, err = quadgk(x -> exp(-x^2), 0, i, rtol=1e-8)
println("Result: $integral")
end
Result: 0.746824132812427
Result: 0.8820813907624216
...
Result: 0.886226925452758
Result: 0.8862269254527579
0.878343 seconds (1.46 M allocations: 87.892 MiB, 2.96% gc time)
```
, so it’s:
a) typical Julia compilation start
b) vscode compilation start
That isn’t a problem for me. Sorry for problem. I will test it more. If it goes about error evaluation can I just avoid assignement of error statement to variable(to not product more unnecessary variables).
---
### Author:  [@Elrod](https://discourse.julialang.org/u/Elrod)
#### Post date: [May 21, 2019, 9:47pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/12 "2019-05-21T21:47:35Z")
> [@wiktorkujawa](#):
>
> @time integral, err = quadgk(x → exp(-x^2), 0, 1, rtol=1e-8)
If you want accurate estimates without the compilation time, you can use `BenchmarkTools`.
```julia
julia> using QuadGK, BenchmarkTools
julia> @time integral, err = quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)
1.194962 seconds (3.32 M allocations: 181.200 MiB, 9.70% gc time)
(0.746824132812427, 7.887024366937112e-13)
julia> @btime quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)
748.218 ns (19 allocations: 592 bytes)
(0.746824132812427, 7.887024366937112e-13)
julia> @benchmark quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)
BenchmarkTools.Trial:
memory estimate: 592 bytes
allocs estimate: 19
--------------
minimum time: 726.023 ns (0.00% GC)
median time: 963.829 ns (0.00% GC)
mean time: 1.014 μs (11.09% GC)
maximum time: 393.248 μs (99.73% GC)
--------------
samples: 10000
evals/sample: 129
```
`@btime` shows the minimum time, while `@benchmark` shows a summary. As you can see, `@benchmark` ran the function 1.29 million times.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 22, 2019, 2:03am UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/13 "2019-05-22T02:03:55Z")
> [@wiktorkujawa](#):
>
> `@time for i=1:100`
[Don’t benchmark in global scope](https://docs.julialang.org/en/latest/manual/performance-tips/#Avoid-global-variables-1) unless you are using `@btime` with [interpolation of globals](https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/doc/manual.md#interpolating-values-into-benchmark-expressions). Also note that integrating from `0` to `i` means that the integration computation (for the same accuracy) gets progressively more expensive as `i` gets larger.
> [@wiktorkujawa](#):
>
> can I just avoid assignement of error statement to variable(to not product more unnecessary variables).
Whether you assign a value to a variable or not makes no difference in performance (especially in a function, where the compiler can just eliminate the variable if it wants to anyway).
---
### Author:  [@wiktorkujawa](https://discourse.julialang.org/u/wiktorkujawa)
#### Post date: [July 13, 2019, 9:41pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/14 "2019-07-13T21:41:42Z")
I have little problem. In Cuda i want to calculate some vector in type of :
```julia
x = (blockIdx().x-1) * blockDim().x + threadIdx().x
y = (blockIdx().y-1) * blockDim().y + threadIdx().y
z = (blockIdx().z-1) * blockDim().z + threadIdx().z
offset=x+y*blockDim().x *gridDim.x+z*blockDim().x *gridDim.x*blockDim().y *gridDim.y
Table[offset]=someValues*quadgk(u->f(u), 0.0f0, upperBound[offset])
```
The problem is that quadgk return a tuple and if I’d want to do it in separate operations like
```julia
integral[offset],error=quadgk(u->f(u), 0.0f0, upperBound[offset])
Table[offset]=someValues*integral[offset]
```
I would need to synchronize threads(or maybe I’m missing something) which isn’t the best idea. Is there an easier way to do it?
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [July 13, 2019, 10:28pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/15 "2019-07-13T22:28:50Z")
> [@wiktorkujawa](#):
>
> The problem is that quadgk return a tuple
You can always do `quadgk(...)[1]` to get the first element of the tuple (the integral) without assigning it to a variable.
Note also that `quadgk(u -> f(u), ...)` is equivalent to `quadgk(f, ...)`
Finally note that, If you are calculating `F(x) = quadgk(f, 0, x)` for a whole bunch of different values of `x ∈ [0,a]`, then there are _much_ more efficient things that you can do than separate `quadgk` integrals for each `x`, which wastes a _lot_ of calculations that could be shared between different `x` values, especially if `f(x)` is a smooth function. For example, using ApproxFun.jl, you can do `F = cumsum(Fun(f, 0..a))`, and then evaluate `F(x)` very quickly — this works by first constructing a polynomial approximation of `f(x)` on `[0,a]` and then forming the polynomial `F(x)` that is the indefinite integral.
---
### Author:  [@wiktorkujawa](https://discourse.julialang.org/u/wiktorkujawa)
#### Post date: [July 28, 2019, 11:28pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/16 "2019-07-28T23:28:55Z")
I have one more problem with calculating inetgral. On CPU I primarily assigned integral=0.0f0 and next used try/catch method to calculate integral, because in some nodes integral can’t be calculated.
The critical part of cuda function looks like that.
```julia
... (
Sx=Segment[XelementIndex+1]-Segment[XelementIndex]; Sy=Segment[YelementIndex+1]-Segment[YelementIndex];
L=sqrt(Sx*Sx+Sy*Sy);
mi=(Segment[ZelementIndex+1]-Segment[ZelementIndex])/Sy;
δz=z[blockIdx().z]+mi*Segment[YelementIndex]-Segment[ZelementIndex];
α=Sx/L; β=Sy/L;
a=1+mi*mi; p=(α*δx+β*δy+mi*δz)/a; q=δx*δx+δy*δy+δz*δz-p*p*a;
integral=0.0f0;
try
integral = quadgk(u -> ((u-p)^2+q)^(-3/2), 0, L, rtol=1f-8)[1]
catch
end;
@atomic Bx[offset]+=integral*(β*δz-mi*δy);
@atomic By[offset]+=integral*(α*δz-mi*δx);
@atomic Bz[offset]+=integral*(α*δy-β*δx))
```
When I delete try/catch integral part, the function works.
With this try/catch quadGK part I got error:
```julia
KernelError: recursion is currently not supported
Try inspecting the generated code with any of the @device_code_... macros.
Stacktrace:
[1] sort! at sort.jl:521 (repeats 2 times)
[2] fpsort! at sort.jl:1106
[3] eignewt at C:\Users\Wiktor\.julia\packages\QuadGK\v6Zrs\src\gausskronrod.jl:43
[4] kronrod at C:\Users\Wiktor\.julia\packages\QuadGK\v6Zrs\src\gausskronrod.jl:137
[5] do_quadgk at C:\Users\Wiktor\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:7
[6] biotSavartCalculation at c:\Users\Wiktor\CudaTests\testfunc.jl:73
```
Is the way to handle with this?
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [August 20, 2019, 12:03pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/17 "2019-08-20T12:03:40Z")
You should notice that your integrand can be evaluated exactly without any quadgk:
\int\_0^L du \frac{1}{\left((u-p)^2+q\right)^{3/2}} = \frac{L-p}{q \sqrt{L^2-2 L p+p^2+q}}+\frac{p}{q \sqrt{p^2+q}}
---
### Author:  [@wiktorkujawa](https://discourse.julialang.org/u/wiktorkujawa)
#### Post date: [August 20, 2019, 2:41pm UTC](https://discourse.julialang.org/t/julia-integral-calculation-community-module-or-own-module/24278/18 "2019-08-20T14:41:31Z")
I know it, it’s really old topic, but thanks for intentions. I’ll close it if it’s possible.