Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: Double Integral error
Replies: 5   Last Post: Nov 28, 2012 1:24 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Francesco Perrone

Posts: 39
Registered: 5/2/12
Double Integral error
Posted: Nov 28, 2012 4:58 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

Good Morning everybody,

I have this code:

clc, clear all, close all
%% Inizialization
time = 614.4; % Analysis Time
Uhub = 11;
HubHt = 90;

TI = 'A'; % Turbulent Intensity (A,B,C as in the IEC or Specific value)
N1 = 512;
N2 = 32;
N3 = 32;

N = N1*N2*N3; % Total Number of Points
t = 0:(time/(N1-1)):time; % Sampled Time Vector

L1 = Uhub*time; % Box length along X
L2 = 150; % Box length along Y
L3 = 150; % Box length along Z

dx = L1/N1; % Grid Resolution along X-axis
dy = L2/N2; % Grid Resolution along Y-axis
dz = L3/N3; % Grid Resolution along Z-axis

V = L1*L2*L3; % Analysis Box Volume

%constants and derived parameters from IEC
gamma = 3.9; %IEC, (B.12)
alpha = 0.2; %IEC, sect. 6.3.1.2
c = 1.476;
b = 5.6;

%set lambda1 according to guidelines (chap.6)
if HubHt < 60
lambda1 = 0.7*HubHt;
else
lambda1 = 42;
end
L = 0.8*lambda1;%IEC, (B.12)

%IEC, Table 1, p.22
if isequal(TI,'A')
Iref = 0.16;
sigma1 = Iref*(0.75*Uhub + b);
elseif isequal(TI,'B')
Iref = 0.14;
sigma1 = Iref*(0.75*Uhub + b);
elseif isequal(TI,'C')
Iref = 0.12;
sigma1 = Iref*(0.75*Uhub + b);
else
sigma1 = str2num(TI)*Uhub/100;
end

%IEC, sect. 6.3.1.3
sigma_iso = 0.55*sigma1;%IEC, (B.12)
sigma2 = 0.7*sigma1;
sigma3 = 0.5*sigma1;


%% Wave number vectors
ik1 = cat(2,(-N1/2:-1/2),(1/2:N1/2));
ik2 = -N2/2:N2/2-1;
ik3 = -N3/2:N3/2-1;

[x y z] = ndgrid(ik1,ik2,ik3);

k1 = reshape((2*pi*L/L1)*x,N1*N2*N3,1);
k2 = reshape((2*pi*L/L2)*y,N1*N2*N3,1);
k3 = reshape((2*pi*L/L3)*z,N1*N2*N3,1);

k = sqrt(k1.^2 + k2.^2 + k3.^2);

%% Calculation of beta by means of the Energy Spectrum Integration
E = @(k) (1.453*k.^4)./((1 + k.^2).^(17/6));
E_integrated = (1.453*k.^4)./((1 + k.^2).^(17/6));
%//Independent integration on segments
Local_int = arrayfun(@(i)quad(E,k(i-1),k(i)), 2:(N1*N2*N3));
%//integral additivity + cumulative removal of queues
E_int = 1.5 - [0 (cumsum((Local_int)))]; %//To remove queues

E_int = reshape(E_int,N,1);
S = k.*sqrt(E_int);
beta = (c*gamma)./S;


%% Help Parameters
k30 = k3 + beta.*k1;

k0 = sqrt(k1.^2 + k2.^2 + k30.^2);

C1 = (beta.*k1.^2.*(k1.^2 + k2.^2 - k3.*k30))./(k.^2.*(k1.^2 + k2.^2));
C2 = (k2.*k0.^2./((k1.^2 + k2.^2).^(3/2))).*atan2((beta.*k1.*sqrt(k1.^2 + k2.^2)),(k0.^2 - k30.*k1.*beta));

xhsi1 = C1 - (k2./k1).*C2;
xhsi2 = (k2./k1).*C1 + C2;

E_k0 = (1.453*k0.^4)./((1 + k0.^2).^(17/6));

C = sigma_iso*sqrt((2*pi^2*L^3/V).*(E_k0./k0.^4));

phi_33 = @(k2,k3) (E_k0./(4.*pi.*k.^4)).*((k1.^2 + k2.^2));

F_33 = arrayfun(@(i) dblquad(phi_33,k3(i),k3(i+1),k2(i),k2(i+1)), 1:((N1*N2*N3)-1));

Unfortunately, when calling F_33, I get the following error msg:

Error using +
Matrix dimensions must agree.

Error in @(k2,k3)(E_k0./(4.*pi.*k.^4)).*((k1.^2+k2.^2))


Error in quad (line 72)
y = f(x, varargin{:});

Error in dblquad>innerintegral (line 82)
Q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:});

Error in quad (line 72)
y = f(x, varargin{:});

Error in dblquad (line 58)
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...

Error in @(i)dblquad(phi_33,k3(i),k3(i+1),k2(i),k2(i+1))

but I don't get the reason of this error, especially cause the dimensions of the involved parameters are consistent.

Could you please shed a light?

I thank you all for supporting.

Best regards,
FPE



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.