This topic may be considered, in some sense, as a follow-up of my previous thread. Thanks to stevengj and ChrisRackauckas I believe I more or less solved the “symbolic part” of my problem and could start working on the actual thing. In short, I need to numerically compute integrals of the kind
f_1(t,t',x) = \int_{t_0}^t \mathrm{d}\tau \int_{t_0}^{t'} \mathrm{d}\tau' \int_0^{\infty} \mathrm{d}y \,F_1(\tau,\tau',y,t,t',x)
and
f_2(t,t',x) = \int_{t}^{t_2} \mathrm{d}\tau \int_{t}^{t'} \mathrm{d}\tau' \int_0^{\infty} \mathrm{d}y \,F_2(\tau,\tau',y,t,t',x),
with analytical expressions for F_1 and F_2 being given, on some grid \lbrace t_i,t_j',x_k \rbrace. For concreteness, let’s say
F_1(\tau,\tau',y,t,t',x) = \cos[x(t - \tau)] \cos[(x + y)(\tau - \tau')]\cos[y(\tau' - t')] (x + y)^{-2} (\tau + t)^{-1} = F_2(\ldots)
As explained in the previous thread, actual expressions are way too lengthy (and also “more singular”, by the way) and are given by some product of other functions and derivatives of thereof, so I use Symbolics.jl to first get the expression and then convert it to a Julia-type function in order to perform the integration. For the latter, I decided to use the HCubature.jl package. So here’s my current attempt of computing f_1 (with some particular choice of t_0 and “infinity”):
using Symbolics
using HCubature
t_grid = collect(0.1:0.1:2);
t_pr_grid = collect(0.1:0.1:2);
x_grid = collect(0.1:0.1:2);
f1 = Array{Float64}(undef,length(t_grid),length(t_pr_grid),length(x_grid));
@variables tau, tau_pr, t, t_pr, x, y;
F_1 = cos(x*(tau - t)) * cos((x + y)*(tau - tau_pr)) * cos(y*(tau_pr - t_pr)) * (x + y)^(-2) * (tau + t)^(-1);
function f1_eval(F_1::Num,t_cur::Real,t_pr_cur::Real,x_cur::Real)
xmin = [1e-6,1e-6,1e-8];
xmax = [t_cur,t_pr_cur,100];
# F_1(tau,tau',y,t,t',x) function at the given grid point (t,t',x):
F1 = substitute(F_1,Dict((t,t_pr,x) .=> (t_cur,t_pr_cur,x_cur)));
F1_func = build_function(F1,[tau,tau_pr,y],expression=Val{false});
return hcubature(F1_func,xmin,xmax,rtol=sqrt(1e-14));
end
for i in 1:length(t_grid)
for j in 1:length(t_pr_grid)
for k in 1:length(x_grid)
f1[i,j,k] = f1_eval(F_1,t_grid[i],t_pr_grid[j],x_grid[k])[1];
end
end
end
Here are my questions:
- This sort of does the job, but I am pretty sure there is a room for improvement. In particular, instead of converting F_1 at each grid point I could maybe convert F_1 to a vector-type function which for a given argument
(tau,tau_pr,y)returns an array F_1(\tau,\tau',y,t_i,t_j',x_k). I am sure that this can be done withSymbolics.jl, but
- I got a bit confused with two types of functions that
build_functionreturns in this case: “…one for out-of-place AbstractArray output and a second which is a mutating function”. I didn’t really understand how to work with either of them, to be honest; - I don’t have fixed boundaries but rather, again, an array of integration limits. Is there any effective way of integrating an array of functions over an array of integration domains? Note that there is a one-to-one correspondence between the integrands and the integration limits.
- So far I don’t have any good idea on what to do with the second integral, where the “temporal integration domain” is no longer a rectangle but a triangle, albeit a particularly simple one. Unfortunately,
HCubature.jldoes not support integration domains of such kind (unless I miss something), but maybe there is some trick to still do this usingHCubature.jl? The obvious one would be to extend the integration domain to a full rectangle and compensate this by modifying F_2 with a bunch of Heaviside functions. But I reckon that’s a very bad idea. Perhaps there is a better one? If not, are there any good Julia packages that support integration over triangular domains?
Thank you in advance!
