Help solving for and/or while loop.

Hi all,

I have a loop which gives me the frequency increments as a final result. These increments are obtained by dividing are under the function into equal area segments. I will copy paste the loop here to make it more clear:

Hs=3;

Tz=8;

Tp=1.408*Tz;

wp=2*pi/Tp;

u=20; % initial number of harmonics w_axis=3;

spacing_w=5000;

w=linspace(0,w_axis,spacing_w);

Area_under_spectrum =@(w) (5/16*Hs^2.*(wp^4./w.^5).*exp(-5/4.*(w./wp).^(-4)));

Area_total=integral(Area_under_spectrum, 0, w_axis);

area=Area_total/u;

for index=2:u+1;

    area_under_frequency_bandwidth=@(w) integral(Area_under_spectrum,0,w);

    g_function=@(w) area_under_frequency_bandwidth(w)-(index-1)*area;

    w(index)=fzero(g_function,[1e-10,w_axis]);

    frequency(index)=w(index);

    frequency_values=frequency(1,1:index);

    frequency_values(1)=0;

    delta_omega(index)=abs(frequency_values(index)-frequency_values(index-1));

end

Now, the problem is that at very low or very high frequencies the frequency increments (delta_omega) is relatively high and I want to limit it to some maximum value. Therefore, at very low or very high frequencies I want the area of a single segment to be obtained in a following manner: if a particular segment requires a frequency increment greater than the maximum I specify I want the segment area to be continuously halved until the required increment is less than the maximum.

At the end I want to have the frequency increments for all segments.

Please help me guys as I cannot solve it ;/

Thanks in advance!

[Geoff Hayes]

Andrejus - in your above code, what is the segment area that you wish to halve? Is that w(index) for the segment given by index, or is that area? And if you do halve some value, what happens to the "unused" area?

[Andrejus]

Hello Geoff,

Thanks for your reaction!

The area segment is area. Well basically what should happen is that there should be two smaller segments and if two are not enough to meet the condition then there should be 4 and so on. Also, I need to know what are the boundary values (frequency) of these segments.

As you can see, the frequency increment at low and high frequency values are large which violate the wave representation. For these locations I want to do that the area would be halved until the frequency increment is below certain value. I think it could be done using for loop and if condition, but I cannot make it... ;/

Answer 1

[Mike Hosea]

In my view this is mostly a software design issue. You probably don't see how to do it because of the way your code is organized in one big thing.

One thing that isn't clear to me is whether you intend to divide the interval [w(i-1),w(i)] evenly to meet the maximum criterion or divide it into a power of 2 segments with equal area. Either can be done easily enough, but the code is different, obviously.

I would approach this with two basic ideas:

1. Take the principles of the code above and write function that computes the next w value with a given w value and the target area for a segment. You probably want this function to be something like wnew = equalAreas(fun,wold,area); In that case the integral call would be based on integral(f,wold,wnew) rather than integral(f,0,wnew), and the function that goes into fzero would subtract the area variable rather than (i - 1)*area.
2. Subdivide the resulting interval if necessary. If the interval is to be evenly divided, you'll need pow2(ceil(log2((w(i+1)-w(i))/maximum))) subintervals (just needs a nested loop). If, rather, the area is to be subdivided evenly, you'll use the equalAreas function with the third input being area/2^m for m = 1, 2, 3, ... until the maximum delta is less than the max.

[Andrejus]

I was trying to implement your 2nd approach but I was not able to... is there a change that you could help me out with the code?

In any case, thank you very much!

[Mike Hosea]

The two ideas I mentioned above are parts of only one approach. Suppose you had this function for finding b such that integral(fun,a,b) == area:

function b = findRightLimit(fun,a,area,searchInterval)
areaFunction = @(b)integral(fun,a,b);
zeroFunction = @(b)areaFunction(fun,a,b) - area;
b = fzero(zeroFunction,searchInterval);

You could call it like this:

w(j) = findRightLimit(Area_under_spectrum,w(j-1),area_total/n,[max(1e-10,w(j-1)),w_axis]

in a while loop that terminates when w(j) >= w_axis. You just start out with w(1) = 0 and in a loop find w(2), w(3), etc.

Now suppose while you are doing that you discover that w(j) - w(j-1) is larger than you like. Well, in that case, as I said, it wasn't clear to me whether you wanted to sub-divide the interval [w(j-1),w(j)] evenly or generate a set of w values that divided the area of the segment equally. If the latter, you would start a nested while loop on k using your findRightLimit function with a target area of area_total/n/2^k and increment k until no w-increment is too large.

There are many minor programming points to work out. For example, if you had to add extra points because an interval w(j) - w(j-1) was too large, you'll have to increment j by 2^k before moving on instead of incrementing j by 1. You could look at it as j = j + 2^k, where k=0 when w(j) - w(j-1) is not too large.

[Andrejus]

It took me some time to realize your proposal but finally it worked. Thank you