MATLAB: Problem with 4-d integral i MatLab

MATLAB

My problem is concerning a 4-d integral. I have read alot around the the forums but as I understand there does not exist a quadquad(@intgn,a1,a2,b1,b2,c1,c2,d1,d2,tol) for a 4-d integral the same way as the tripelquad(@intgn,a1,a2,b1,b2,c1,c2,tol) in MatLab. I have seen Peter Spellucci respons which is posted more than 1 place around the forums concerning this problem, but I cant make it work. My function is:
function Int = intgn(r,z1,z2,z) 
lm=0.00625;
s=0.00625;
%z=1e-3;
k_lms=sqrt(4*(lm+s)*r.*(((lm+s)+r).^2+(z-z1)^(-1)));
k_s=sqrt(4*(s)*r.*(((s)+r).^2+(z-z1)^(-1)));
K_lms=mfun('EllipticK',sqrt(k_lms));
K_s=mfun('EllipticK',sqrt(k_s));
E_lms=mfun('EllipticE',sqrt(k_lms));
E_s=mfun('EllipticE',sqrt(k_s));
konst=1;
F_lms1=(z-z1)./(r.*((lm+s+r).^2+(z-z1)^2).^(1/2)).*(-K_lms+((lm+s)^2+r.^2+(z-z1).^2).*E_lms./((lm+s-r).^2+(z-z1)^2));
F_s1=(z-z1)./(r.*((s+r).^2+(z-z1)^2).^(1/2)).*(-K_s+((s)^2+r.^2+(z-z1)^2).*E_s./((s-r).^2+(z-z1)^2));
F1=2*(F_lms1-F_s1);
F_lms2=(z-z2)./(r.*((lm+s+r).^2+(z-z2)^2).^(1/2)).*(-K_lms+((lm+s)^2+r.^2+(z-z2)^2).*E_lms./((lm+s-r).^2+(z-z2)^2));
F_s2=(z-z2)./(r.*((s+r).^2+(z-z2)^2).^(1/2)).*(-K_s+((s)^2+r.^2+(z-z2)^2).*E_s./((s-r).^2+(z-z2)^2));
F2=2*(F_lms2-F_s2);
Int=konst.*r.*F1.*F2;
———————————————————— And the main file is: ————————————————————
tic
clear all 
clc
r_in=6.25e-3;
r_out=12.5e-3;
tau_m=5e-3;
%z=tau_m/2
%Magnet
mu_0=pi*4e-7;
M0=891000 %Fra Bachelor
sigma=3.82e7; %Conductors conduction (alu)
I=M0*tau_m;
tol=1e4;
tripvaerdi = triplequad(@intgn, r_in, r_out, -tau_m, tau_m, -tau_m, tau_m, tol);
toc
The tripelint works fine as long as I remove the '%' before z=1e-3, but I cant make the file integrate over all 4 variable.
Does anyone got any suggestions?
Best regards
Jonas

Best Answer

You can insert a single integral as an argument in triplequad:
zmin=-tau_m; zmax=tau_m;
f1 = @(r,z1,z2) quadv(@(z) intgn(r,z1,z2,z),zmin,zmax,tol);
quadInt = triplequad(f1, r_in, r_out, -tau_m, tau_m, -tau_m, tau_m, tol);
It is important to use quadv for the single integral because you need a vector output for each scalar input.
However, you should keep in mind the caution in this thread.
Related Question