Nested differentiation in nested integral

Hello everyone!

I have quite the complex integral that looks like this:

I am trying to evaluate this function, but I'm having issues.

Where this gets troublesome, is the d/dr differentation nested in between an integral2 and integral3 loop.

Here is my code so far:

n = 179;
ID = 1e-3; %c - w/2 in figure
OD = 2e-3; %c + w/2 in figure
h = 1e-3;
R = 1e-3;
d = 1e-3;
Z = 1e-3;
c = 1e-3;
I = 2;
J = 1;
u0 = 1;
M = 1;

Fz = funtest(u0, J, M, R, d, ID, OD, Z, h)

function [Fz] = funtest(u0, J, M, R, d, ID, OD, Z, h)
    syms r
    f1 = @(theta, rho, phi, z, r) rho./((z - d/2).^2 + r.^2 + rho.^2 - 2.*r.*rho.*cos(theta - phi)).^(1/2);
    f2 = @(theta, rho, phi, z, r) rho./((z + d/2).^2 + r.^2 + rho.^2 - 2.*r.*rho.*cos(theta - phi)).^(1/2);
    
    f1Q = @(phi, z, r) integral2(@(theta, rho) f1(theta, rho, phi, z, r), 0, 2 * pi, 0, R, 'AbsTol', 1e-5, 'RelTol', 1e-3);
    f2Q = @(phi, z, r) integral2(@(theta, rho) f2(theta, rho, phi, z, r), 0, 2 * pi, 0, R);

    fQ = @(phi, z, r) (f1Q(phi, z, r) - f2Q(phi, z, r));
    fQdr = @(phi, z, r) diff(fQ(phi, z, r), r) * r; % Line causing trouble, how to best tackle?

    f = integral3(@(phi, z, r) arrayfun(fQdr, phi, z, r), 0, 2 * pi, Z - h/2, Z + h/2, ID/2, OD/2, 'AbsTol', 1e-5, 'RelTol', 1e-3);
    Fz = J*u0*M*f/(4*pi);
end

fQdr = @(phi, z, r) diff(fQ(phi, z, r), r) * r;

Without the above line, the function seems to work well (i.e. I use fQ instead of fQdr in the integral3 function) and commenting the fQdr line out.

I am clueless how I should tackle this complex integral including the differentation.

Does anyone have any idea how I could tackle this? It'd be super appreciated.

Thank you!

Update:

I've managed to achieve something, but it's not right yet. I now perform the integration and differentiation manually. The code is slow and also gives from output, any idea how to properly implement the equation?

function [Fz] = funtest(u0, J, M, R, d, ID, OD, Z, h)
    
    r = linspace(ID, OD, 10);
    dr = r(2) - r(1);

    fr = zeros(length(r), 1);
    for i = 2:length(r)
        f1 = @(theta, rho, phi, z) rho./((z - d/2).^2 + r(i).^2 + rho.^2 - 2.*r(i).*rho.*cos(theta - phi)).^(1/2);
        f2 = @(theta, rho, phi, z) rho./((z + d/2).^2 + r(i).^2 + rho.^2 - 2.*r(i).*rho.*cos(theta - phi)).^(1/2);

        f1_1 = @(theta, rho, phi, z) rho./((z - d/2).^2 + r(i - 1).^2 + rho.^2 - 2.*r(i - 1).*rho.*cos(theta - phi)).^(1/2);
        f2_2 = @(theta, rho, phi, z) rho./((z + d/2).^2 + r(i - 1).^2 + rho.^2 - 2.*r(i - 1).*rho.*cos(theta - phi)).^(1/2);
        
        f1Q = @(phi, z) integral2(@(theta, rho) f1(theta, rho, phi, z), 0, 2 * pi, 0, R, 'AbsTol', 1e-5, 'RelTol', 1e-3);
        f2Q = @(phi, z) integral2(@(theta, rho) f2(theta, rho, phi, z), 0, 2 * pi, 0, R, 'AbsTol', 1e-5, 'RelTol', 1e-3);

        f1Q_1 = @(phi, z) integral2(@(theta, rho) f1_1(theta, rho, phi, z), 0, 2 * pi, 0, R, 'AbsTol', 1e-5, 'RelTol', 1e-3);
        f2Q_2 = @(phi, z) integral2(@(theta, rho) f2_2(theta, rho, phi, z), 0, 2 * pi, 0, R, 'AbsTol', 1e-5, 'RelTol', 1e-3);
    
        fQ = @(phi, z) ((f1Q(phi, z) - f2Q(phi, z)) - (f1Q_1(phi, z) - f2Q_2(phi, z)))./dr.*r(i);
    
        fr(i) = integral2(@(phi, z) arrayfun(fQ, phi, z), 0, 2 * pi, Z - h/2, Z + h/2, 'AbsTol', 1e-5, 'RelTol', 1e-3);
    end

    frtable = ones(length(r) - 1, 1);
    for i = 1:(length(r) - 1)
        frtable(i) = (fr(i + 1) - fr(i))/2 .* dr;
    end
    f = sum(frtable(1:end-1));
    Fz = J*u0*M*f/(4*pi);

end

Would appreciate the help a ton!

Thanks!

Answer 1

[Torsten]

I would (without mathematical scruple) differentiate with respect to "r" under both integrals, then use "integral2" to integrate both (differentiated) functions with respect to "theta" and "rho" and pass the sum of the results of these integrations to "integral3".

[Timothy]

Hello kind stranger!

Thorsten, as your title names you, you're the true MVP.

Thank you, I didn't realize I could just move the differentiator. It works like a charm and I get the output I require.

Thanks a ton!!

[David Goodmanson]

Hi Timothy,

One thing you can do right off before you do anything else, is reduce the 5 dimensional integral to a 4 dimensional one. The theta integration is nested inside the phi integration, therefore phi is a constant for the theta integration. You can let theta = sigma + phi, the integrand now involves cos(sigma) and the limits of the sigma integration are -phi to 2pi-phi. Since you are going all the way around the circle with a trig function you can call the limits 0 to 2pi instead. You end up with an integral dsigma from 0 to 2pi of an integrand involving cos(sigma). Now nothing in the entire integrand involves phi any more, so the dphi integration just gives a factor of 2pi, leaving a 4 dimesional integral.

[Timothy]

Hi David!

Thank you for this insight! I'm definitely going to try this and see if I can improve on computational speed! Very much appreciated!

Best, Timothy