Numerical integration methods for surface integral

40 次查看（过去 30 天）

显示更早的评论

Lorenzo Lellini 2023-2-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1909045-numerical-integration-methods-for-surface-integral

在 MATLAB Online 中打开

Hello. I have to perform this integral:

and this other one that is basically the same:

where:

a is a known complex vector that is the FFT of the derivative of a real laboratory mesuerement
are known
S is a circumference surface with radius R (known)

Since a is a complex numerical vector, I tried to use cumtrapz and integral functions to compute this integral. The two methods give me worng results and the results are different depending on the used method.

Moreover, I have seen that the final result of the method 2 with the integral function is the same for the first integral (pa_a) and for the second one (p).

Could you please help me to find the right way to perform this calculation?

%% DATA
% Input
y = data;
% Compute acceleration from measured velocity
a = diff(y);
A = fft(a);
% Measurements data
Fs = 44100;     % Measure sampling frquency
f_start = 2;    % Initial measure frequency [Hz]
f_end = 5000;   % Final measure frequency [Hz]
% Parameters
rho_0 = 1.225; % kg/m3
ra = 1; % m 
rc = 0; % m 
p_0 = 1; % Standard pressure
c = 340; % Sound speed m/s
R = 0.254;   % Radius
Sc = pi*R^2; % Membrane surface [m2]
f = linspace(f_start, f_end, length(diff(fft(y)))); % Frequency vector
k = 2*pi*f/c; % Wavenumber
%% Integral computation
% METHOD 1: trapezoidal integration
domain1 = linspace(0, 2*pi, length(A));
domain2 = linspace(0, R, length(A));
pa_a_cumtrapz = rho_0/(2*pi) * cumtrapz(domain1, cumtrapz(domain2, A/abs(ra-rc)));
p_cumtrapz = rho_0/(2*pi) * cumtrapz(domain1, cumtrapz(domain2, A/abs(ra-rc) .* exp(-1i*k*abs(ra-rc))));
%------------------------
% METHOD 2: Double integral (I cannot use integral2 because I have a vector and not a matrix)
% First internal integral computation
integrand = @(om) A/abs(ra-rc);
temp = integral(integrand, 0, R, "ArrayValued",true);
% Second integral computation
integrand2 = @(om) temp;
pa_a_integral = rho_0/(2*pi).* integral(integrand2, 0, 2*pi, "ArrayValued",true);
p_integral = rho_0/(2*pi) * integral(integrand2, 0, 2*pi, "ArrayValued",true) .* exp(-1i*k*abs(ra-rc));
%% INTEGRATION RESULT PLOTS
figure;
loglog(abs(pa_a_cumtrapz));
hold on
loglog(abs(p_cumtrapz)*0.3) % *0.3 is a shift to check differences between two spectra
legend ('pa_a cumtrapz', 'p cumtrapz low shifted')
xlim([f_start f_end]);
grid minor
%------------------------
figure;
loglog(abs(pa_a_integral));
hold on
loglog(abs(p_integral)*0.3) % *0.3 is a shift to check differences between two spectra
legend ('pa_a integral', 'p integral low shifted')
xlim([f_start f_end]);
grid minor