# how to speed up integration

10 visualizaciones (últimos 30 días)
Dominika el 2 de Mayo de 2014
Comentada: Dominika el 5 de Mayo de 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 comentariosMostrar NingunoOcultar Ninguno
Roger Stafford el 5 de Mayo de 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.
Dominika el 5 de Mayo de 2014
it's 78 degrees indeed! Thanks for this remark Roger! My problem is solved now :)

Iniciar sesión para comentar.

Friedrich el 2 de Mayo de 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 comentariosMostrar 8 comentarios más antiguosOcultar 8 comentarios más antiguos
Friedrich el 5 de Mayo de 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.
Dominika el 5 de Mayo de 2014
OK, Thanks a lot!

Iniciar sesión para comentar.

### Categorías

Más información sobre Programming en Help Center y File Exchange.

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by