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!
采纳的回答
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".
3 个评论
Timothy
2023-5-25
David Goodmanson
2023-5-25
编辑：David Goodmanson
2023-5-25
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
2023-5-26
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Translated by 
