how to speed up integration
7 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Hi,
I do integration in my code:
% fun= @(teta) ((abs((((-1i.*(T.*(k0.*sin(teta)).^2+EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2))./(omega.*Zair))+((1i.*(EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2))./(omega.*Zair)))./((1-((1i.*(EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2))./(omega.*Zair))).*(1+((-1i.*(T.*(k0.*sin(teta)).^2+EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2))./(omega.*Zair)))))).^2).*sin(teta).*cos(teta)
tau_tr_dash= integral(fun,0,78,'ArrayValued',true);
The process takes more than 10min!
How can I speed it up?
My whole code:
function sandwich_unidirectional (rho_c,d,rho,EI,c,b,vc,Ec)
f = [31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000....
1250 1600 2000 2500 3150 4000 5000 6300 8000]; % frequency range [Hz]
Zair=428 ; % air impedance at 0°C [N*s/m^3]
%rho_c-core material density
M0=(rho_c*b*c)/6;
E0=(Ec*10^9)/(1-vc^2);
v0=vc/(1-vc);
K=(E0*b)/(c*(1-v0^2));
A=b*d; %-cross secional area of the face sheet
%d-thickness of the face sheet
%rho-density of sheet plate material
%EI-flexural rigidity of the face sheet
omega=2.*pi.*f; % angular frequency [rad/s]
c_air=331.6; % sound speed in air at 0° [m/s]
k0=omega./c_air;
%k=k0.*sin(teta)
%c-thickness of the core
%b-hight of the panel
%vc-Poisson's ratio of the core
%Ec-elastic modulus of the core
T=(E0*b)/(c*(1-v0^2));
% fun= @(teta) ((abs((((-1i.*(T.*(k0.*sin(teta)).^2+EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2))./(omega.*Zair))+((1i.*(EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2))./(omega.*Zair)))./((1-((1i.*(EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2))./(omega.*Zair))).*(1+((-1i.*(T.*(k0.*sin(teta)).^2+EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2))./(omega.*Zair)))))).^2).*sin(teta).*cos(teta)
tau_tr_dash= integral(fun,0,78,'ArrayValued',true);
TL=10.*log10(abs(1./tau_tr_dash))
semilogx(f,TL,'k');
grid on
ylabel('Transmission Loss [dB]');
xlabel ('Frequency')
title('Transmission Loss')
hold on
end
Thanks for any suggestions,
Dominika
2 Commenti
Roger Stafford
il 5 Mag 2014
Are you integrating from 0 to 78 radians, or is it from 0 to 78 degrees? The way you have it set up is for 0 to 78 radians and that would take you through about twenty-five full cycles in sin(teta) which would require quite a bit more processing than if it were only 78 degrees.
Risposta accettata
Friedrich
il 2 Mag 2014
Hi,
Can you provide values for all your variables? S what is EI,k0 etc.
What normally helps to speed integration up is to use MATLAB Coder to generate a MEX file out of the integration part. Integral itself is not supported for code generation but quadgk is. If you can provide some values for your variable I could test it and post the timing result.
What do you set 'ArrayValued' to true? As far as I see its not a vector valued function. fun seems like a R^1 => C^1 function.
Have you tested other integral functions like quadgk?
10 Commenti
Friedrich
il 5 Mag 2014
Non of the integration routines give a correct result. All give you a wrong answer because of the nature of your fun. It's not easy to integrate that numerically because of the singularities fun has. You would either need to adjust fun to make it a smooth function or check some litature of numerical integration to see if there any methods which can deal with that kind of functions.
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!