Numerical integration of vector-valued functions using Julia packages

Specific Domains Numerics
1

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

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!

2

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 provides a function expv(t, A,b) that computes exp(t*A)*b.

For actual timing,

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)
3 Likes
3

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: I think you’re interested in computing H(Δ).

2 Likes
4

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):

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

gives

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> 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.)

13 Likes