how to speed up integration

Question

Dominika 2014-5-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/127962-how-to-speed-up-integration

评论： Dominika 2014-5-5

采纳的回答： Friedrich

在 MATLAB Online 中打开

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.

Dominika 2014-5-5

it's 78 degrees indeed! Thanks for this remark Roger! My problem is solved now :)

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Friedrich 2014-5-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/127962-how-to-speed-up-integration#answer_135365

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 个评论
显示 8更早的评论隐藏 8更早的评论

Friedrich 2014-5-2

在 MATLAB Online 中打开

Do you know the result you would expect? Because if I integrat that fun with using quadv i get a completly different result and when I split up that arrayvalued function into scalar valued functions and integrate them seperatly I also get a different result.

When using your code (so integral function)

 TL =
Columns 1 through 9
   31.1860   33.0053   34.7157   36.5106   38.4013   40.2040   42.0404   44.1052   45.9946
Columns 10 through 18
   47.8998   49.8846   51.9444   53.8734   55.8746   57.9455   59.8812   61.8177   63.9607
Columns 19 through 25
   65.8982   67.8359   69.8430   71.9178   73.8558   75.8631   77.9380

When replacing integral with quadv I get (way faster)

 TL =
Columns 1 through 9
   57.3525   66.5634   75.4774   84.9571   94.9493  104.4064  113.9439  124.5568  134.1862
Columns 10 through 18
143.8375  153.8484  164.2063  173.8873  183.9178  194.2884  203.9769  213.6663  224.3862
Columns 19 through 25
234.0766  243.7672  253.8040  264.1788  273.8697  283.9066  294.2815

And when splitting it up by using this code (also way faster than integral with array valued set to true)

function TL = sandwich_unidirectional()
%inputs hardcoded
rho_c=20;
d=0.001;
rho=7800;
EI=210;
c=0.06;
b=2;
vc=0.3;
Ec=0.03;
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));
TL = zeros(numel(f),1);
for j=1:numel(f)
     omega=2.*pi.*f(j);
     k0=omega./c_air;
     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 = quadgk(fun,0,78);
     tau_tr_dash = integral(fun,0,78);
     TL(j)=10.*log10(abs(1./tau_tr_dash));
end
% 
% semilogx(f,TL,'k');
% grid on
% ylabel('Transmission Loss [dB]');
% xlabel ('Frequency')
% title('Transmission Loss')
% hold on
end

I get

At least when replacing integral with quadgk the result stays the same. Can you shed some light onthe result you would expect?

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.

Dominika 2014-5-5

OK, Thanks a lot!

请先登录，再进行评论。

how to speed up integration

2 个评论
显示无隐藏无

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

how to speed up integration

2 个评论 显示 无隐藏 无

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

10 个评论
显示 8更早的评论隐藏 8更早的评论