# how to speed up integration

26 views (last 30 days)
Dominika on 2 May 2014
Commented: Dominika on 5 May 2014
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
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 CommentsShowHide 1 older comment
Dominika on 5 May 2014
it's 78 degrees indeed! Thanks for this remark Roger! My problem is solved now :)

Friedrich on 2 May 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?
Dominika on 5 May 2014
OK, Thanks a lot!