how to speed up integration

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

[Roger Stafford]

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]

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

Answer 1

[Friedrich]

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]

Thanks for an answer Friedrich,

The integral is a vector 1x25 and that's why I set ArrayValued to true. I tried different ways but there is always problem with 'matrix dimentions must agree'. I added the whole code to my main post.

Thanks a lot!

[Friedrich]

I am still missing the inputs to your sandwich_unidirectional function. What is rho_c,d,rho,EI,c,b,vc,Ec?

[Dominika]

aaa yeah.. so sorry..

rho_c=20
d=0.001
rho=7800
EI=210
c=0.06
b=2
vc=0.3
Ec=0.03

Thank you :)

[Friedrich]

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

TL = 
   31.1860
   33.0053
   34.7157
   36.5106
   38.4013
   40.2040
   42.0404
   44.1052
   45.9946
   47.8998
   49.8846
   93.1474
   90.7580
   73.1180
   57.9455
   59.8812
  104.3994
  119.7487
  118.1098
  107.1495
  120.9006
  108.1523
  123.8067
  136.8159
  148.8706

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

[Friedrich]

I used some symbolic math on it to get a better impression of fun and it seems like fun is pretty evil when it comes to get inegrated because it has singularities and scales horribly. So the overal numerical integration you do on that fun is basically the same as simply calling rand(25,1) to get a result. I dont see anything you can do here to get a reliable result.

[Dominika]

results obtained using my code seem to be OK to me. But results from quadv and splitting the code are way too high.

[Dominika]

So what do you suggest me to do in this case? How can I use 'rand' in my case?

[Dominika]

and actually the 'fun' expression is a little bit changed, it should be:

fun= @(teta) ((abs((((-1i.*(T.*(k0.*sin(teta)).^2+EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2+2.*K-M0.*omega.^2))./(omega.*Zair))+((1i.*(EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2+3.*T.*(k0.*sin(teta)).^2-3.*M0.*omega.^2))./(omega.*Zair)))./((1-((1i.*(EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2+3.*T.*(k0.*sin(teta)).^2-3.*M0.*omega.^2))./(omega.*Zair))).*(1+((-1i.*(T.*(k0.*sin(teta)).^2+EI.*10.^9.*(k0.*sin(teta)).^4-rho.*A.*omega.^2+2.*K-M0.*omega.^2))./(omega.*Zair)))))).^2)

the results are slightly changed, but probably it doesn't affect the solution of the problem

[Friedrich]

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]

OK, Thanks a lot!