How do you create a loop to find a coefficient value to which one of the variables in the function changes with every cycle

For example:

The change in weight (dW) of an aircraft due to fuel consumption is shown below:

       dW=(-SFC/N)*sqrt((2*W^3)/(p*S))*(1/((Cl^1.5)/Cd))*dt

With 'W' being the total weight of the aircraft and current fuel load, 'dt' being the change in time of the aircraft and the rest of the coefficients being constant. For a given time period of 100,000 seconds, how would the loop be generated? As the 'W' variable would change with every cycle, resulting in a different dW in every cycle ad so on. I believe the following should be involved but I'm not sure how, any help would be greatly appreciated.

           dt=0:100:100000
           W=W+dW
           dW=(-SFC/N)*sqrt((2*W^3)/(p*S))*(1/((Cl^1.5)/Cd))*dt

Thanks alot James

[Oleg Komarov]

You don't need a loop. Plug in all the values...and whenever you need an element by element operation between SAME size vectors, use the .() operators (ex: .*, ./).

Or post all the relevant values for the vars involved and the full error message if you don't succeed in the implementation.

Answer 1

[David Young]

The first thing to be clear about is the distinction between the time step, dt, and the time itself, which might be called t. dt is the change in time between successive values of t.

You probably want dt to be constant, but t will change. So your loop might look something like this:

t_start = 0;
t_end = 100000;
dt = 100;  % time step

for t = t_start : dt : t_end
  % compute dW and update W here
end

[James Norris]

Thanks David, I've corrected my distinction of the time step, but my main difficulty is with computing the values of dW with respect to the newly obtained values of W from every cycle?

which are then to be plotted later

Thanks for helping

[David Young]

James, I think you just use the formulae that you showed in your original post to compute dW and W inside the loop.

Answer 2

[Walter Roberson]

The code shown is independent of the time if the time-step is constant. The code can then be rewritten as a differential equation. In Maple notation, it would be:

dsolve({W(0) = W0, diff(W(t), t) = -SFC*sqrt(2*W(t)^3/(p*S))*Cd/(N*Cl^1.5)})

where W0 is the initial weight.

Maple returns the following solution:

W(t) = RootOf(2*t*SFC*Cd+N*sqrt(2)*Cl^(3/2)*(Int(1/sqrt(_a^3/(p*S)), _a = _b .. _Z))-N*sqrt(2)*Cl^(3/2)*(Int(1/sqrt(_a^3/(p*S)), _a = _b .. W0)))

I must admit that I have no idea what the introduced variable _b is supposed to mean. If I ask Maple to combine the parts, I get

W(t) = RootOf(Int(N*sqrt(2)*Cl^(3/2)/sqrt(_a^3/(p*S)), _a = W0 .. _Z)+2*t*SFC*Cd)

which basically means that W(t) is the value of the variable _Z that causes the expression inside the RootOf() to evaluate to 0. This will be a decreasing value as t increases.

The Int() part can be evaluated to a closed form under some circumstances, depending upon the signs of the values involved. If the known constraints on the constants allow for it, then a ratio of two polynomials results, and equating to -2*t*SFC*Cd and solving for _Z gives an exact ratio-of-polynomials form and thus an exact form for W(t) without needing simulation.

[David Young]

An instructive exercise is to compare the results from the Euler method numerical integration with the results from the analytical solution, and particularly observing how changes to delta-t affect the errors.