Using nested for loops to achieve double summation

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

Answer 1

[Sand Man]

Anyone has an idea?