Effacer les filtres
Effacer les filtres

Integrate function by for loop

4 vues (au cours des 30 derniers jours)
Kaiyuan
Kaiyuan le 4 Mai 2024
Commenté : Kaiyuan le 6 Mai 2024
clear
Cs1=4; Cs2=5.6; Cl1=2; Cl2=2.8; ns1=1.5; ns2=-1.5; nl1=1.5; nl2=-1.5;
rho=1.22; A=10; V=44.72; l=20; n=1;Ua=20;
yaw=0:0.5:90;
t=1:1:60;
u=[4.16886273000000 -5.46160601200000 -2.21536928300000 -0.107926889000000 4.17995704900000 2.70294649300000 -3.57644741500000 2.00451722300000 0.881703288000000 4.62541932200000 1.63637915200000 5.73084754400000 0.961007086000000 0.350184264000000 1.37768570200000 -7.64386640500000 1.09384013000000 -3.32210770200000 1.06977546700000 3.77699477800000 -2.80761673400000 -2.51814521100000 -2.56997163700000 -0.434123529000000 -7.33456835300000 3.18164556400000 -2.46259887500000 -2.81705263700000 -6.11824449800000 4.20959983600000 0.958003023000000 1.34911158700000 -3.66536927800000 1.33702769200000 0.951579011000000 4.81211100200000 1.70504507900000 7.08168644200000 3.20839901100000 -0.438521707000000 2.92055007900000 -0.280320940000000 -1.25512163200000 -2.68273256600000 -5.07343734300000 -1.42260189600000 4.56907055100000 3.27623428200000 -1.23514645600000 -2.86948142600000 1.75847520500000 2.79193783800000 -4.18180552400000 -3.14073219500000 -2.43496557500000 -1.28117719800000 4.95079755900000 -3.52809836900000 1.96375591300000 -5.66592848200000]
for i=1:length(yaw)
ms(i)=2*sin(yaw(i)*pi/180)*V/l;
ml(i)=2.5*sin(yaw(i)*pi/180)*V/l;
if yaw(i)<=40
Cs(i)=Cs1*(sin(yaw(i)*pi/180)/sin(30*pi/180))^ns1;
Cl(i)=Cl1*(sin(yaw(i)*pi/180)/sin(30*pi/180))^nl1;
elseif yaw(i)>40
Cs(i)=Cs2*(sin(yaw(i)*pi/180)/sin(30*pi/180))^ns2;
Cl(i)=Cl2*(sin(yaw(i)*pi/180)/sin(30*pi/180))^nl2;
end
S1(i)=0.5*rho*A*(V^2)*Cs(i); % Steady side forces (Uniform along each row)
L1(i)=0.5*rho*A*(V^2)*Cl(i); % Steady lift forces (Uniform along each row)
Msi=ms(i);
Mli=ml(i);
Csi=Cs(i);
Cli=Cl(i);
Yawi=yaw(i*pi/180);
for j=1:length(t)
Tj=t(j);
Uj=u(j);
fun1=@(t1)(((2.*pi.*n.*Msi)^2).*t1.*exp(-2.*pi.*n.*Msi.*t1)).*Uj.*(Tj-t1);
fun2=@(t1)(((2.*pi.*n.*Mli)^2).*t1.*exp(-2.*pi.*n.*Mli.*t1)).*Uj.*(Tj-t1);
S2(i,j)=rho*A*Csi*Ua*(1+(1/(2*Csi))*diff(Csi)*cot(Yawi))*integral(fun1,0,inf); % Unsteady side forces (noniform along each row and column)
L2(i,j)=rho*A*Csi*Ua*(1+(1/(2*Csi))*diff(Csi)*cot(Yawi))*integral(fun2,0,inf); % Unsteady lift forces (noniform along each row and column)
S(i,j)=Sl(i)+S2(i,j); % Total side forces
L(i,j)=S1(i)+S2(i,j); % Total lift forces
j=j+1
end
i=i+1;
end
Dear experts:
The equations I am going to solve are attached as above screenshot and the code is listed.
To summary, the equation is the: Total force(i,j)=Mean force (i)+Unsteady force (i,j (involve intergration over time interval t1)). Hence I will get i*j number of total force result. It is 181*60 (row*column) data.
The mean force componant S1, L1 is uniform value inside each row, which mean that the mean force does not vary along columnn on each row (i); but the mean force componants only has different value between rows.
The calculation of the unsteady force S2,L2 is more complex and the unsteady force is different between each row and column (i,j), and inside each (i,j), the unsteady force is calculated based on the intergration over time interal (0,infinite):S2 (i,j)=(eq...)*integration(fun1(t1), 0 ,inf), L2 (i,j)=(eq...)*integration(fun2(t1), 0 ,inf),
Thereforce, the total force S(i,j)=S1(i)+S2(i,j), L(i,j)=L1(i)+L2(i,j)
The code does not have problem with calculating the mean force (i) S1, L1, but it shows "Array indices must be positive integers or logical values" after I add the for loop (j) to calculate the unsteay force (i,j) S2, L2 which involve integration over time interval t1.
Could our expert help me fix the code?
Kind Regards

Réponses (1)

John D'Errico
John D'Errico le 4 Mai 2024
Modifié(e) : John D'Errico le 4 Mai 2024
As always look carefully at your code. Learn to use the debugger. It would tell you the problem, and in which line there was a problem. I scanned through your code, and I saw the error, just waiting to be found. But we can do it the hard way too.
Cs1=4; Cs2=5.6; Cl1=2; Cl2=2.8; ns1=1.5; ns2=-1.5; nl1=1.5; nl2=-1.5;
rho=1.22; A=10; V=44.72; l=20; n=1;Ua=20;
yaw=0:0.5:90;
t=1:1:60;
u=[4.16886273000000 -5.46160601200000 -2.21536928300000 -0.107926889000000 4.17995704900000 2.70294649300000 -3.57644741500000 2.00451722300000 0.881703288000000 4.62541932200000 1.63637915200000 5.73084754400000 0.961007086000000 0.350184264000000 1.37768570200000 -7.64386640500000 1.09384013000000 -3.32210770200000 1.06977546700000 3.77699477800000 -2.80761673400000 -2.51814521100000 -2.56997163700000 -0.434123529000000 -7.33456835300000 3.18164556400000 -2.46259887500000 -2.81705263700000 -6.11824449800000 4.20959983600000 0.958003023000000 1.34911158700000 -3.66536927800000 1.33702769200000 0.951579011000000 4.81211100200000 1.70504507900000 7.08168644200000 3.20839901100000 -0.438521707000000 2.92055007900000 -0.280320940000000 -1.25512163200000 -2.68273256600000 -5.07343734300000 -1.42260189600000 4.56907055100000 3.27623428200000 -1.23514645600000 -2.86948142600000 1.75847520500000 2.79193783800000 -4.18180552400000 -3.14073219500000 -2.43496557500000 -1.28117719800000 4.95079755900000 -3.52809836900000 1.96375591300000 -5.66592848200000]
u = 1x60
4.1689 -5.4616 -2.2154 -0.1079 4.1800 2.7029 -3.5764 2.0045 0.8817 4.6254 1.6364 5.7308 0.9610 0.3502 1.3777 -7.6439 1.0938 -3.3221 1.0698 3.7770 -2.8076 -2.5181 -2.5700 -0.4341 -7.3346 3.1816 -2.4626 -2.8171 -6.1182 4.2096
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
for i=1:length(yaw)
ms(i)=2*sin(yaw(i)*pi/180)*V/l;
ml(i)=2.5*sin(yaw(i)*pi/180)*V/l;
if yaw(i)<=40
Cs(i)=Cs1*(sin(yaw(i)*pi/180)/sin(30*pi/180))^ns1;
Cl(i)=Cl1*(sin(yaw(i)*pi/180)/sin(30*pi/180))^nl1;
elseif yaw(i)>40
Cs(i)=Cs2*(sin(yaw(i)*pi/180)/sin(30*pi/180))^ns2;
Cl(i)=Cl2*(sin(yaw(i)*pi/180)/sin(30*pi/180))^nl2;
end
S1(i)=0.5*rho*A*(V^2)*Cs(i); % Steady side forces (Uniform along each row)
L1(i)=0.5*rho*A*(V^2)*Cl(i); % Steady lift forces (Uniform along each row)
Msi=ms(i);
Mli=ml(i);
Csi=Cs(i);
Cli=Cl(i);
Yawi=yaw(i*pi/180);
for j=1:length(t)
Tj=t(j);
Uj=u(j);
fun1=@(t1)(((2.*pi.*n.*Msi)^2).*t1.*exp(-2.*pi.*n.*Msi.*t1)).*Uj.*(Tj-t1);
fun2=@(t1)(((2.*pi.*n.*Mli)^2).*t1.*exp(-2.*pi.*n.*Mli.*t1)).*Uj.*(Tj-t1);
S2(i,j)=rho*A*Csi*Ua*(1+(1/(2*Csi))*diff(Csi)*cot(Yawi))*integral(fun1,0,inf); % Unsteady side forces (noniform along each row and column)
L2(i,j)=rho*A*Csi*Ua*(1+(1/(2*Csi))*diff(Csi)*cot(Yawi))*integral(fun2,0,inf); % Unsteady lift forces (noniform along each row and column)
S(i,j)=Sl(i)+S2(i,j); % Total side forces
L(i,j)=S1(i)+S2(i,j); % Total lift forces
j=j+1
end
i=i+1;
end
Array indices must be positive integers or logical values.
Now, where are you indexing an array or vector? I found it pretty quickly. What is yaw? yaw is a vector. What do you do with it?
Yawi=yaw(i*pi/180);
Is i*pi/180 a positive integer? I think that is unlikely, given that pi is a transcendental number. What you wanted to do there, I don't know. Maybe you wanted to index the vector yaw to find the ith element, then convert to radians? I cannot guess. Your code, not mine. I'll gues, however, that you WANTD to write:
Yawi = yaw(i)*pi/180;
Since you do something similar in other places.
There may be other errors of this sort, but that was the only one that jumps out.
  1 commentaire
Kaiyuan
Kaiyuan le 6 Mai 2024
Hi John:
I appreciate your patience to read the code and reply.
I will follow the advice and try to identify the errors.
Regards

Connectez-vous pour commenter.

Catégories

En savoir plus sur Loops and Conditional Statements dans Help Center et File Exchange

Produits


Version

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by