IF loop within a FOR loop problem

[kieranunited9]

I am editing code which was developed to model a mass spring damper system in the time domain using the equations of motion of the system. One damping ratio in particular,cpto, was constant originally but I am trying to vary it according to the if statement below.

if k3*(x2(i)-x1(i))>=0 && k3*(x2(i)-x1(i))<=rel
cpto(i+1)=inf;
elseif k3*(x2(i)-x1(i))<=0
cpto(i)=0;
else cpto(i)=2.1*sqrt(k3*mc);

end


The if statement is inside a larger for loop which determines the output of the system over 10000 iterations. The problem is that x2 and xc are used to determine cpto. For this reason I want to either use past values of x2(i-1) and xc(i-1) to determine a current value of cpto(i) (at the start of the for loop) or else use current values of x2(i) and xc(i) to determine a future value of cpto(i+1) at the end of the loop.

The if statement works if put at the end of the for loop, if it calculates cpto for current values of x2 and xc but the problem is it's too late for the loop to use cpto.

Any ideas why it won't work when I use future values of cpto or why i can't use past values of x2 and xc, they are in use elsewhere in the for loop anyway. Thanks for your time,

Kieran.

 

[vivek.iitd]

It would be helpful if you can be more clear. what do you mean by
"The if statement works if put at the end of the for loop, if it calculates cpto for current values of x2 and xc but the problem is it's too late for the loop to use cpto."
and where is xc ?

 

[kieranunited9]

Thanks for the reply.
Apologies, I wasn't very clear. I have posted a version of the for loop below. All variables are initialised and initial conditions specified. I'm not sure why it's not working, I have decided to use the past values of x2 and xc to determine cpto for the current loop. The output is crazy though. I am a novice at this, how do I go about debugging it? Would a while loop be dificult to implement instead of the if statement?

for i = 2:n

if (k3*(x2(i-1)-xc(i-1))>=0) && (k3*(x2(i-1)-xc(i-1))<=rel)
cpto(i)=999999999999999;

elseif k3*(x2(i)-xc(i))<=0
cpto(i)=0;

else cpto(i)=2.1*sqrt(k3*mc);

endl(i) = amp*sin(w1*t(i-1));

F(i) = (k1*(x0(i-1)-x1(i-1)))+(c1*(v0(i-1)- v1(i-1)));

Fint(i) = Fint(i-1) + F(i)*dt;

r(i) = 0.5*(x0(i-1)-(1/z1)*Fint(i));

x0(i) = l(i)+ r(i); % Equations of motiona1(i) = ((-c1*(v1(i-1)-v0(i-1)))-(k1*(x1(i-1)-x0(i-1)))-(k2*(x1(i-1)-x2(i-1)))-(c2*(v1(i-1)-v2(i-1))))/m1;

v1(i) = v1(i-1) + a1(i)*dt;

x1(i) = x1(i-1) + v1(i)*dt;

a2(i) =((-c2*(v2(i-1)-v1(i-1)))-(k2*(x2(i-1) - x1(i-1)))-(c3*(v2(i-1)-vc(i-1)))-(k3*(x2(i-1)-xc(i-1)))-(kl*(x2(i-1)-x5(i-1))))/m2;

v2(i) = v2(i-1) + a2(i)*dt;

x2(i) = x2(i-1) + v2(i)*dt;

ac(i)= ((-c3*(vc(i-1)-v2(i-1)))-(c4*(vc(i-1)-v4(i-1)))-(k3*(xc(i-1)-x2(i-1)))-(k4*(xc(i-1)-x4(i-1)))-((cpto(i))*vc(i-1)))/mc;

vc(i) = vc(i-1) + ac(i)*dt;

xc(i) = xc(i-1) + vc(i)*dt;

ForceWEC(i)= ((-c3*(vc(i-1)-v2(i-1)))-(c4*(vc(i-1)-v4(i-1)))-(k3*(xc(i-1)-x2(i-1)))-(k4*(xc(i-1)-x4(i-1))));a4(i) = ((-c4*(v4(i-1)-vc(i-1)))-(k4*(x4(i-1)-xc(i-1)))-(cshore*v4(i-1))-(kl*(x4(i-1)-x5(i-1))))/m4;

v4(i) = v4(i-1) + a4(i)*dt;

x4(i) = x4(i-1) + v4(i)*dt;% Leakage

a5(i) = ((-kl*(x5(i-1)-x2(i-1)))-(kl*(x5(i-1)-x4(i-1))))/m5;

v5(i) = v5(i-1) + a5(i)*dt;

x5(i) = x5(i-1) + v5(i)*dt;


v0(i) = (x0(i)-x0(i-1))/dt;
vl0(i) = (l(i)-l(i-1))/dt;
vr0(i) = (r(i)-r(i-1))/dt;
t(i) = t(i-1) + dt;
end

 

[Unrest]

You want to update cpto each iteration, and use that updated value in the next iteration?

Why is cpto an array? Why not a single value that changes every iteration? Do you need to keep track of all the previous values?

Either way, you can update cpto at the end of the loop and the new value you put in will be ready for the next iteration. Doing it that way means you have to initialize cpto with something so the first iteration has some value to start with.

eg:

cpto(1) = initial value
for i = 1 to n
 x2(i)=1234*cpto(i)
 cpto(i+1) = 4321*x2(i)
end for

In your original code you seem to be mixing up cpto(i) and cpto(i+1) in the same "if" statement.

Incidentally, it's bad practice to use variables names like "i". Call it "timeStep" or whatever it actually means, so nobody will have to struggle to interpret it. Same goes for cpto, k3, x2, etc. This kind of extreme shorthand is convenient when you're writing equations on paper, but unnecessary and confusing in code.

 

[vivek.iitd]

I see a problem in the "If" statement, it uses past values of x2 and xc i.e. x2(i-1) and xc(i-1)
as well as present values of x2 and xc i.e. x2(i) and xc(i) for the current loop i.


------------------------------------------------------------------------------
if (k3*(x2(i-1)-xc(i-1))>=0) && (k3*(x2(i-1)-xc(i-1))<=rel)
cpto(i)=999999999999999;

elseif k3*(x2(i)-xc(i))<=0
cpto(i)=0;

else cpto(i)=2.1*sqrt(k3*mc);

end
------------------------------------------------------------------------------

now there are two cases:
1."If" coming at the beginning.
2."If" coming at the end.

As you said earlier that "If" works fine at the end of the loop, but not at the beginning.
There is a thing to notice here, that you are also changing the current values of x2 and xc i.e. x2(i) and xc(i) according to x2(i-1) and xc(i-1).

------------------------------------------------------------------------------
x2(i) = x2(i-1) + v2(i)*dt;
xc(i) = xc(i-1) + vc(i)*dt;
------------------------------------------------------------------------------


Therefore, if you use "If" statement at the end, then it will use NEW current values of x2 and xc [modified as above] and past value of x2 and xc.


Where as in the other case when the "If" statement comes at the beginning, it will use OLD current values of x2 and xc (as there are modified later in the code). Therefore, in this case the answer will be completely different from the case when "If" is used at the end.

Hope this helps.


Also, i agree with Unrest, that if you are not required to keep track of all the values of cpto, then it will be better to define it as a single variable changing values for each iteration.

 

[vivek.iitd]

Or, may be "If" statement should look like this



if (k3*(x2(i-1)-xc(i-1))>=0) && (k3*(x2(i-1)-xc(i-1))<=rel)
cpto(i)=999999999999999;

elseif k3*(x2(i-1)-xc(i-1))<=0
cpto(i)=0;

else cpto(i)=2.1*sqrt(k3*mc);

end