Problem with 4-d integral i 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

Answer 1

[Andrew Newell]

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.

[Jonas]

Thank you very much - It really helped.