I am trying to perform a numerical triple integral over s, gamma1, gamma2. The limits are (-inf +int), (0,+inf) and (gamma1,+inf) respectively. Please dont be scared from the shape of my function (its just a function of gamma1, gamma2, s)
The following is my code
syms s
syms gamma1
syms gamma2
fun=-(exp(-(28035689158432973*pi*gamma2^(2/3))/2305843009213693952)*
exp(-(pi*s*7120816246010697*i)/112589990684262400)*
(1/((pi*s*(4194304/gamma1^2 + 4194304/gamma2^2)*i)/(50*
(6144/gamma1 + 6144/gamma2)) + 1)^((3*(2048/gamma1 + 2048/gamma2)^2)
/(4194304/gamma1^2 + 4194304/gamma2^2)) - 1)
*(exp(-(pi^2*s*(log((-(gamma2*25*i)/(1024*pi*s))^(1/3) + 1)/3;
y=@(s,gamma1,gamma2)fun;
gamma2min=@(s,gamma1) gamma1;
prob=integral3(y,-inf,+inf,0,+inf,gamma2min,+inf)
I get the following error
Error using integralCalc/finalInputChecks (line 511)
Input function must return 'double' or 'single' values. Found 'sym'.
Any advice?
Thank you very much!
You can use quadgk to numerically integrate functions defined like
y=@(a,b,c) 1/abs(a^2+b^2+c^2+1);
(which I used to test my answer).
It's tricky because quadgk expects a function that takes vector input and returns a vector of function values, but you can get around it by using a lot of arrayfun's:
R=@(s,gamma1) quadgk(@(gamma2) arrayfun(@(k) y(s,gamma1,k),gamma2),gamma1,Inf)
S=@(s) quadgk(@(gamma1) arrayfun(@(k) R(s,k),gamma1),0,Inf)
T=quadgk(@(s) arrayfun(@(k) S(k),s),-Inf,Inf)
But! It's very slow, I wasn't patient enough to wait for the answer. So, replace the Inf and -Inf limits with 100 and -100, for example, and you will get an answer. Maybe try with 50 and -50 and see how much the solution changes by, if it changes very little then you can be confident the answer is quite accurate, otherwise you'll have to increase the number and wait longer! The faster your function decays the smaller the bounds you will be ale to get away with.