how to speed up integration
7 次查看(过去 30 天)
显示 更早的评论
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 个评论
Roger Stafford
2014-5-5
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.
采纳的回答
Friedrich
2014-5-2
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 个评论
Friedrich
2014-5-5
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.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!