Help solving for and/or while loop.

1 view (last 30 days)
Andrejus on 4 Nov 2014
Commented: Andrejus on 13 Nov 2014
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.
Andrejus on 4 Nov 2014
Hello Geoff,
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... ;/

Mike Hosea on 5 Nov 2014
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 on 13 Nov 2014
It took me some time to realize your proposal but finally it worked. Thank you