Numerical integration methods for surface integral
36 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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
0 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differentiation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!