# Numerical integration of vector-valued functions using Julia packages

**URL:** <https://discourse.julialang.org/t/numerical-integration-of-vector-valued-functions-using-julia-packages/85413>\
**Category:** Numerics\
**Tags:** integral\
**Created:** [August 6, 2022, 7:59pm UTC](https://discourse.julialang.org/t/numerical-integration-of-vector-valued-functions-using-julia-packages/85413 "2022-08-06T19:59:26Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![danrib07](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danrib07/32/212751_2.png) [@danrib07](https://discourse.julialang.org/u/danrib07)\
**Post date:** [August 6, 2022, 7:59pm UTC](https://discourse.julialang.org/t/numerical-integration-of-vector-valued-functions-using-julia-packages/85413/1 "2022-08-06T19:59:26Z")

</div>

Hi,

I need to compute the following integral \int\_{0}^{t} e^{At'}b \,dt' where A is a real-valued matrix of size (m x m), and b is a real-valued column vector of size (m x 1). The current method I have implemented to compute this integral is using the package `Trapz` such as to compute the integral for every row of e^{At'}b. Here’s a working example

```julia
using Trapz

function f(A,b,t)
    return exp(A * t)*b
end

t_init = 0
t_final = 10
step = 1
time = t_init:step:t_final

A = [1 2; 3 4]
b = [1;2]

integrand = zeros(2,length(time))
for (idx,t) in enumerate(time)
    integrand[:,idx] = f(A,b, t)
end

result = zeros(2,1)
for n = 1:2
    result[n] = trapz(time,integrand[n,:])
end

```

This seems to work, but I was wondering if there is a more concise/faster way to compute this integral using one of the integration packages currently available in Julia? Any helps is greatly appreciated it!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 6, 2022, 8:29pm UTC](https://discourse.julialang.org/t/numerical-integration-of-vector-valued-functions-using-julia-packages/85413/2 "2022-08-06T20:29:50Z")

</div>

The biggest thing you can do is speed up `f`. It’s possible to use iterative methods to compute this without ever storing `exp(A*t)` which are much faster for larger `A` and `b`. For example, [https://github.com/SciML/ExponentialUtilities.jl](https://ExponentialUtilities.jl) provides a function `expv(t, A,b)` that computes `exp(t*A)*b`.

For actual timing,

```julia
A=rand(1000,1000);
b=rand(1000);
@time exp(3.0 *A)*b;
  0.175631 seconds (25 allocations: 53.781 MiB)
@time expv(3.0, A,b);
  0.004479 seconds (27 allocations: 282.844 KiB)

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [August 6, 2022, 9:06pm UTC](https://discourse.julialang.org/t/numerical-integration-of-vector-valued-functions-using-julia-packages/85413/3 "2022-08-06T21:06:57Z")

</div>

If I remember correctly, these integrals can also be computed by taking a block from the exponential of a larger matrix. Take a look at this [paper by Van Loan](https://www.olemartin.no/artikler/vanloan.pdf): I think you’re interested in computing `H(Δ)`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 6, 2022, 10:26pm UTC](https://discourse.julialang.org/t/numerical-integration-of-vector-valued-functions-using-julia-packages/85413/4 "2022-08-06T22:26:32Z")

</div>

> [@danrib07](#):
>
> I need to compute the following integral \int\_{0}^{t} e^{At'}b \,dt' […] The current method I have implemented to compute this integral is using the package `Trapz`

First, I would almost _never_ use the trapezoidal rule to compute numerical integrals of smooth integrands when you can evaluate the integrand wherever you want. It is _exponentially_ faster to use something like Gaussian quadrature. Moreover, packages like QuadGK can handle vector-valued integrands directly (rather than integrating each row separately):

```julia
using QuadGK
A = [1 2; 3 4]
b = [1;2]
t_final = 10
quadgk(t -> exp(A*t)*b, 0, t_final)[1]

```

gives

```julia
2-element Vector{Float64}:
 3.7347743661398145e22
 8.164742103848492e22

```

which is accurate to at least 12 digits.

However, in this particular case, the integral is _exactly_ A^{-1} (e^{At} - I) b:

```julia
julia> using LinearAlgebra

julia> A \ ((exp(A*t_final) - I) * b)
2-element Vector{Float64}:
 3.7347743661399164e22
 8.164742103848739e22

```

(Of course, you can speed this up further as noted above by @Oscar_Smith by using a package to compute e^{At} b directly rather than computing e^{At} first.)