Using nested for loops to achieve double summation

2 views (last 30 days)
Hey there,
I am trying to calculate the available power in irregular waves in the time domain. The code I have written is based on the following formula:
To give more insight into the formula, the first single summation in the formula is understood followingly. At certain time t all the regular wave components from n=1 to N are added together to form the irregular wave. And this is repeated for the chosen time length. The same principle applies for the other part in the formula, but instead it is a double summation.
There is a way to check if the code returns the correct result. If the simulation is done for a long enough ( e.g. 300s) duration, then the average of the calculated power should be the same as calculated by the time average power of the same formula, which is:
I am certain, that the calculation result for the time average power (second formula) is correct. But the the average of the first formula is much higher - up to a power of 10 higher. I am quite certain that I have a mistake somewhere within the nested loops. But I haven't been able to figure it out. Can someone point out where I have made a mistake?
The result for the first formula was: average_timedomain = 4.076230132823717e+04 And the for the second formula: ave = 4.850139959777149e+03
I have also attached the one subfunction I use in my code.
And here is the code itself:
%%Inputs %%
Ohm = [0.01:0.01:4];
Hs = 2.5;
Tp = 9.6;
Gamma = 3.3;
h = 100;
t=[0:1:200];
rho_w = 1025;
g = 9.81;
k = linspace(0.0017,1.631,400);
%%Values from other functions%%
% Jonswap spectrum values %
[S, Amp, Phase]=JONSWAP(Ohm, Hs, Tp);
%%Calculation of the power in irregular waves as a function of time %%
sum=0 %Initialization
for ii = 1:length(Ohm)
sum = sum+rho_w*g* (Amp(ii)*Ohm(ii)/k(ii)^2) * ((cosh(k(ii)*h)+1)/sinh(k(ii)*h)) * ...
cos(k(ii)-Ohm(ii)*t+Phase(ii))
for jj = 1:length(Ohm)
if jj == ii
sum = sum + rho_w*g*((Amp(ii)*Amp(jj)*Ohm(jj))/(cosh(k(ii)*h)*sinh(k(jj)*h))) * ...
(h/2 + ((sinh(h*(k(ii)+k(jj))))/(2*(k(ii)+k(jj)))) ) * ...
cos(k(ii)-Ohm(ii)*t+Phase(ii)) .* cos(k(jj)-Ohm(jj)*t+Phase(jj))
else
sum = sum + rho_w*g*((Amp(ii)*Amp(jj)*Ohm(jj))/(cosh(k(ii)*h)*sinh(k(jj)*h))) * ...
( ((sinh(h*(k(ii)-k(jj))))/(2*(k(ii)-k(jj)))) + ((sinh(h*(k(ii)+k(jj))))/(2* (k(ii)+k(jj)))) ) * ...
cos(k(ii)-Ohm(ii)*t+Phase(ii)) .* cos(k(jj)-Ohm(jj)*t+Phase(jj))
end
end
end
absolut=abs(sum);
average_timedomain=mean(absolut);
%%Calculation of the time average power in irregular waves %%
ave = 0; %Initialization
for nn = 1:length(Ohm)
ave = ave + ((0.25 * rho_w*g*Amp(nn)^2*Ohm(nn))/k(nn)) * (1 + ((2*k(nn)*h)/(sinh(2*k(nn)*h))))
end

Answers (1)

Sand Man
Sand Man on 20 Jan 2014
Anyone has an idea?

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!