I am trying to evaluate a triple integral using the 'trapz' command in MATLAB. My integrand is function of gamma1,gamma2,s... I want to evaluate the integral numerically and not symbolically.
syms gamma1; syms gamma2; syms s
f=gamma1*gamma2*exp(-s*10);
x= -20:0.3:20;
y=linspace(0.1,100000,length(x));
z=linspace(0.1,100000,length(x));
[s,gamma1,gamma2] = ndgrid(x,y,z);
mat=eval(f).* (gamma2>gamma1); %MY QUESTION IS HERE
out = trapz(x,trapz(y,trapz(z,mat,3),2),1);
My question is I have a condition on my integrand, it should be evaluated ONLY when gamma2 >gamma1. Is my code correct above i.e is the way I add the logical statement correct?
Yes, the logic is correct: mat will have zeros where the condition gamma2>gamma1 fails.
You can test correctness yourself by using an example where you know the answer already. For instance, integration of f=gamma1*gamma2*s within [0,2] in each variable gives 8 (with paper and pencil), and if the condition gamma2>gamma1 is added, the result is halved (by symmetry), so it becomes 4. And indeed, your code returns approximately 4 for this function.
As an aside: you may want to reconsider the desirability of using the same number of sample points in every dimension, given the disparity of sizes (40 by 100000 by 100000).