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!

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.

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

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]

However, in this particular case, the integral is exactlyA^{-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.)