Numerical integration methods for surface integral
26 次查看(过去 30 天)
显示 更早的评论
Hello. I have to perform this integral:
and this other one that is basically the same:
where:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1289440/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1289445/image.png)
- 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
0 个评论
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!