Function handle not rechecking if statement

2 Ansichten (letzte 30 Tage)
Ryan O'Connor
Ryan O'Connor am 23 Apr. 2019
Kommentiert: dpb am 28 Apr. 2019
I am trying to model a buoy in sinesoidal waves. I am using function handles to do this. There is 4 if conditions:
  1. When there is no itersection between the buoy and the wave and the buoy is obove water
  2. When there is no itersection between the buoy and the wave and the buoy is below water
  3. When there is itersection between the buoy and the wave less that 50%
  4. When there is no itersection between the buoy and the wave greater than 50%
However what occurs is that it runs through the inital values and determines which if statement to use in the for loop and then just uses that value, even though the values effecting the if statement above the for loop are changing after each loop.
tf_f = @(T,Y) isempty(p_f(T,Y)); % Is the matrix empty == true(1) or false (0) [ is there intersection points from the buoy and wave surface]
%% Conditions for forces
% float
if tf_f(T,Y) == 1 && bottom_f(Y) > p21top_f(T) %if intersection matrix is empty and bottom of buoy is above the wave surface
elseif tf_f(T,Y) == 1 && top_f(Y) < p21top_f(T) %if intersection matrix is empty and top of buoy is below the wave surface
elseif tf_f(T,Y) == 0 && centre_check_f(Y) > p21top_f(T) %if intersection matrix is full and centre of buoy is above the wave surface
elseif tf_f(T,Y) == 0 && centre_check_f(Y) < p21top_f(T) %if intersection matrix is full and centre of buoy is below the wave surface
end
Above is the if statement that decides the how it should act, tf_f is used to determine whether the intersection matrix is empty, however it only does the above for the inital and does not update? after this the Runga_kutta 4th order intergration takesplace in a for loop
%Float (centre)
Accel_h_f = @(T,Y,V) ((Sum_forces_h_f(T,Y,V))/(abs(mass_f)));
for i=1:(length(t)-1) % calculation loop
k_1vh_f = Accel_h_f(t(i), y_f(i), v_f(i))*h;
From what i can tell, all other parts of my code are correct and for each of the conditions above. I just cant seem to figure out how to get the code to check the if statement each time. Any suggestions?

Antworten (1)

dpb
dpb am 23 Apr. 2019
Bearbeitet: dpb am 24 Apr. 2019
tf_f = @(T,Y) isempty(p_f(T,Y));
if tf_f(T,Y) == 1 && bottom_f(Y) > p21top_f(T)
elseif tf_f(T,Y) == 1 && top_f(Y) < p21top_f(T)
elseif tf_f(T,Y) == 0 && centre_check_f(Y) > p21top_f(T)
elseif tf_f(T,Y) == 0 && centre_check_f(Y) < p21top_f(T)
end
Not enough code to be able to tell...we don't know where the loop is which is crucial...
But, the crystal ball is back from the shop so let's give it a whirl and see if it is (finally!) repaired or not--
It says, since the content of any variable not passed as a dummy argument into an anonymous function is the value for that variable in the workspace at the time the anonymous function is defined and is invariant, that definition is outside the loop and so the array p_f is constant in the function definition instead of being revised as you expect.
Either pass p_f, too, or don't use the anonymous function but just write the test expression.
(Now to wait and see if I wasted money on the repair bill... :) )
  2 Kommentare
Ryan O'Connor
Ryan O'Connor am 28 Apr. 2019
Bearbeitet: dpb am 28 Apr. 2019
Apologies for the late repsonse; I was busy with exams. I tried moving the for loop above the if statement and it worked as far as it now changes between the conditions.
However, there is now a seperate issue. While the buoy and the wave are intersecting tf_f = 0 the interx function doesnt appear to find these points causing it to error out. It seems that it fails on the last step in before transitioning, when it is about to transittion from tf=0 to tf = 1.
Index in position 2 exceeds array bounds.
Error in Floating_perfect_Save_Spare>@(Q,Z,S)Q(Z,S)
Error in Floating_perfect_Save_Spare>@(T,Y)select(p(T,Y),1,2)
Although if you look at the plot you can clearly see that there is still intersection so I can't seem to determine why its missing finding these point??
I've tried including a try and catch function. Which to my knowlegde is if it errors do this but it never did the action.
Any idea what could be causing the issue or anyways around it? I've included a good portion of the code so you can see it for yourself since I left out a lot last time
%% Preparation
clc
clear all
close all
%% inital Conditions to be used in the function handles
h=0.005; % step size
t = 0:h:3; % Calculates upto y(3)
y = zeros(1,length(t));
v = zeros(1,length(t));
y(1) = 0;
v(1) = 0;
Y = y(1);
T = t(1);
v = v(1);
%% Constants
g = 9.81; % gravity [m/s2]
DSw = 1025; % Density sea water [kg/m3]
L = 1; % wavelength [m]
A = 0.0375; % wave amplitude [m]
Ti = 1; % Wave Period [s]
c = L/Ti; % wave speed (celerity)[m/s]
d = 0.45; % Water Depth [m]
H = 2*A; %wave height [m]
k = (2*pi)/L; % Wave Numebr
w = (2*pi)/Ti; % Wave Frequency [rads/s]
wl = 0; % water line for wave velocity and acceleration
S_k = 90; % Spring constant [N/m]
%% Coefficents
% Sphere
Cd = 0.1; % Drag
eh = 0.5; %radiation damping heave
es = 0.63; %radiation damping surge
amh = 0.55; % added mass heave
ams = 0.63; % added mass surge
%% Buoy parameters
r=0.075; % radius of buoy
Vbuoy_Total = 4/3*pi*r^3; % total volume of the buoy assuming spherical [m3]
Denpercent = 0.0635;
Den_buoy = DSw*Denpercent;
mass = Den_buoy*Vbuoy_Total; % mass of buoy[kg]
Weight = mass*g; % weight of buoy [N]
Fbmax = DSw*g*Vbuoy_Total; % max possible bupoyancy fully submerged [N]
numwave = 2; % number of desired waves on plot
inc = 150; % Number of increment for step size for time and wavelength
inct = Ti/inc;
ts = 0:inct:numwave*Ti; % time [s]
inclambda = L/inc; % increment step size for a single wave length
ts = 0:inclambda:numwave*L;
th = 0:pi/((numwave*inc)/2):2*pi; % phase angle
%% for loop begin
for i=1:(length(t)-1) % calculation loop
%% Buoy centre and equation
cx = 1; % centre of gravity poisition in x-direction
cy = 0.035 + d; % centre of gravity position in y-direction
xc = r * cos(th) + cx ; % x values of buoy relative to poition in cartesian co-ordinates
yc = @(Y) r * sin(th) + cy + Y; % y values of buoy relative to poition in cartesian co-ordinates
centre_check = @(Y) cy + Y;
crc = @(Y) [xc;yc(Y)]; % Buoy geometry
top = @(Y) cy + r + Y; % highest point of buoy
bottom = @(Y) cy - r + Y; % lowest point of buoy
%% Wave surface elevation
xp = 0:inclambda:numwave*L; % horizontal distance [m] ** make a rounder number**
N =@(T) d + A*cos(k*(xp+cx) + w*T); % Surface of the wave (assuming Airy's linear wave theory)
wprofile = @(T) [xp ; N(T)]; % wave x and y co-ordinate
% water particle velocity
U1_s = @(T) ((pi*H)/Ti)*exp((2*pi*wl)/L)*cos(w*T); % water particle velocity in the surge (x-direction)
U2_h = @(T) ((pi*H)/Ti)*exp((2*pi*wl)/L)*sin(w*T); % water particle velocity in the heave (y-direction)
%% intersection of buoy and wave
%line
x_value(1:4) = cx;
dy = 0.2;
y_value = 0.2:dy:0.8;
line = [x_value;y_value];
select = @(Q,Z,S) Q(Z,S);
ptop = @(T) InterX(line,wprofile(T));
p11top = @(T) select(ptop(T), 1, 1);
p21top = @(T) select(ptop(T), 2, 1);
p =@(T,Y) InterX(crc(Y),wprofile(T)); % intersection points of wave and buoy done using the interX function
select = @(Q,Z,S) Q(Z,S);
tf = @(T,Y) isempty(p(T,Y)); % Is the matrix empty == true(1) or false (0) [ is there intersection points from the buoy and wave surface)
if tf(t(i),y(i)) == 1 && bottom(y(i)) > p21top(t(i)) %if intersection matrix is empty and bottom of buoy is above the wave surface
Fbuoy_h = @(T,Y)(T*Y)*0; % not in wave therefore no buoyancy
F_drag_h = @(T,Y,V) (T*Y*V)*0;
F_drag_h_wp = @(T,Y) (T*Y)*0;
Vol_sub = @(T,Y) (T*Y)*0;
F_Rd_h = @(T,Y)(T*Y)*0;
F_s_h = @(Y) S_k*(centre_check(Y)); % spring force in heave direction
Vol_percentage = @(T,Y) (T*Y)*0;
ifx(i) = 1
elseif tf(t(i),y(i)) == 1 && top(y(i)) < p21top(t(i))
A_sub =@(T,Y) pi*(r^2);
Vol_sub = @(T,Y) Vbuoy_Total;
Fbuoy_h = @(T,Y) (T*Y)*0 + Fbmax;
F_Rd_h = @(T,Y) (1/2)*Vol_sub(T,Y)*w*eh;
F_s_h = @(Y) S_k*(centre_check(Y)); % spring force in heave direction
F_drag_h = @(T,Y,V) 1/2*Cd*DSw*((V)*abs(V))*A_sub(T,Y);
F_drag_h_wp = @(T,Y) 1/2*Cd*DSw*(U2_h(T)*abs(U2_h(T)))*A_sub(T,Y);
Vol_percentage = @(T,Y) (T*Y)*0 +100;
ifx(i) = 2
elseif tf(t(i),y(i)) == 0 && centre_check(y(i)) > p21top(t(i)) %if intersection matrix is full and centre of buoy is above the wave surface
p11 = @(T,Y) select(p(T,Y), 1, 1);
p12 = @(T,Y) select(p(T,Y), 1, 2);
p21 = @(T,Y) select(p(T,Y), 2, 1);
p22 = @(T,Y) select(p(T,Y), 2, 2);
Distance_interX = @(T,Y) sqrt(((p12(T,Y) - p11(T,Y))^2) + ((p22(T,Y) - p21(T,Y))^2)); % distance between intersection point, Spherical cap diameter (2r)
Spherical_cap = @(T,Y)Distance_interX(T,Y)/2; % spherical cap base radius
Spherical_Distance =@(T,Y) sqrt((r^2) - (Spherical_cap(T,Y)^2)); % distance from centre of buoy to wave surface [m]
Cap_height =@(T,Y) r - Spherical_Distance(T,Y); % Spherical cap height relative to wave surface [m]
Vol_sub = @(T,Y)((pi*(Cap_height(T,Y)^2))/3)*((3*r)-Cap_height(T,Y)); % volume of spherical buoy submerged [m3]
Vol_percentage = @(T,Y)Vol_sub(T,Y)/Vbuoy_Total *100; % percentage of buoy submerged
A_sub =@(T,Y)pi*(Spherical_cap(T,Y)^2);
Fbuoy_h = @(T,Y) DSw*g*Vol_sub(T,Y); % buoyancy force expericenced by float [N]
F_drag_h = @(T,Y,V) 1/2*Cd*DSw*((V)*abs(V))*A_sub(T,Y);
F_Rd_h = @(T,Y) (1/2)*Vol_sub(T,Y)*w*eh;
F_s_h = @(Y) S_k*(centre_check(Y)); % spring force in heave direction
F_drag_h_wp = @(T,Y) 1/2*Cd*DSw*(U2_h(T)*abs(U2_h(T)))*A_sub(T,Y);
ifx(i) = 3
elseif tf(t(i),y(i)) == 0 && centre_check(y(i)) < p21top(t(i)) %if intersection matrix is full and centre of buoy is below the wave surface
p11 = @(T,Y) select(p(T,Y), 1, 1);
p12 = @(T,Y) select(p(T,Y), 1, 2);
p21 = @(T,Y) select(p(T,Y), 2, 1);
p22 = @(T,Y) select(p(T,Y), 2, 2);
Distance_interX = @(T,Y) sqrt(((p12(T,Y) - p11(T,Y))^2) + ((p22(T,Y) - p21(T,Y))^2)); % distance between intersection point, Spherical cap diameter (2r)
Spherical_cap = @(T,Y)Distance_interX(T,Y)/2; % spherical cap base radius
Spherical_Distance =@(T,Y) sqrt((r^2) - (Spherical_cap(T,Y)^2)); % distance from centre of buoy to wave surface [m]
Cap_height =@(T,Y) r - Spherical_Distance(T,Y); % Spherical cap height relative to wave surface [m]
Vol_sub = @(T,Y) ((Vbuoy_Total/2) + ((Vbuoy_Total/2) - ((pi*(Cap_height(T,Y)^2))/3)*((3*r)-Cap_height(T,Y)))); % volume of spherical buoy submerged [m3]
A_sub =@(T,Y)pi*(Spherical_cap(T,Y)^2);
Vol_percentage = @(T,Y)Vol_sub(T,Y)/Vbuoy_Total *100; % percentage of buoy submerged
Fbuoy_h = @(T,Y) DSw*g*Vol_sub(T,Y); % buoyancy force expericenced by float [N]
F_drag_h = @(T,Y,V) 1/2*Cd*DSw*((V)*abs(V))*A_sub(T,Y);
F_Rd_h = @(T,Y) (1/2)*Vol_sub(T,Y)*w*eh;
F_s_h = @(Y) S_k*(centre_check(Y)); % spring force in heave direction
F_drag_h_wp = @(T,Y) 1/2*Cd*DSw*(U2_h(T)*abs(U2_h(T)))*A_sub(T,Y);
ifx(i) = 4
end
added_mass_h = @(T,Y) DSw*Vol_sub(T,Y); %added mass for acceleration
Sum_forces_h = @(T,Y,V) Fbuoy_h(T,Y) - Weight + F_drag_h(T,Y,V);% + F_Rd_h(T,Y) + F_s_h(Y) + F_drag_h_wp(T,Y); % + added_mass_y(i);
%% Equation of motion to be intergrated using the Runge-Kutta 4th order method for Ordinary Diffferential Equation's (ODE's)
Accel_h = @(T,Y,V) ((Sum_forces_h(T,Y,V))/(abs(mass + added_mass_h(T,Y))));
k_1vh = Accel_h(t(i), y(i), v(i))*h;
k_1h = (v(i))*h;
k_2vh = Accel_h(t(i) + 0.5*h, y(i) + 0.5*k_1h, v(i) + 0.5*k_1vh)*h;
k_2h = (v(i) + (k_1vh/2))*h;
k_3vh = Accel_h((t(i) + 0.5*h),(y(i) + 0.5*k_2h) , v(i) + 0.5*k_2vh)*h;
k_3h = ((v(i)+(k_2vh/2))*h);
k_4vh = Accel_h((t(i)+h), (y(i)+k_3h), (v(i)+k_3vh))*h;
k_4h = (v(i) + k_3vh)*h;
n(i) = (1/6)*(k_1h + 2*k_2h + 2*k_3h + k_4h);
q(i) = (1/6)*(k_1vh + 2*k_2vh + 2*k_3vh + k_4vh);
y(i+1) = y(i) + n(i); % main equation cx
v(i+1) = v(i) + q(i); % main equation cy
%% Plots
vols(i) = Vol_percentage(t(i),y(i));
fsum(i) = Sum_forces_h(t(i),y(i),v(i));
tfx (i) = tf(t(i),y(i));
p21topp(i) = p21top(t(i));
bottomp(i) = bottom(y(i));
ip11(i) = p11(t(i), y(i));
ip12(i) = p12(t(i), y(i));
ip21(i) = p21(t(i), y(i));
ip22(i) = p22(t(i), y(i));
N = d + A*cos(k*(xp) + w*t(i)); % Surface of the wave (assuming Airy's linear wave theory)
xcn = r * cos(th) + cx; % x values of buoy relative to poition in cartesian co-ordinates
ycx =yc(y(i)); % y values of buoy relative to poition in cartesian co-ordinates
crcx = [xcn;ycx]; % Buoy geometry
figure(1)
title('motion of buoy on wave')
plot(xp,N,...
crcx(1,:),crcx(2,:))
ylim([0.25;0.7])
end
dpb
dpb am 28 Apr. 2019
Well, with all the code you did show, you didn't show us anything related to the error itself...namely the function that had the error.
However since the crystal ball does seem to have been repaired, and it even has a little bit of information to work with this time; let's put it to the test again-- :)
Index in position 2 exceeds array bounds.
Error in Floating_perfect_Save_Spare>@(Q,Z,S)Q(Z,S)
Error in Floating_perfect_Save_Spare>@(T,Y)select(p(T,Y),1,2)
Since the error message says you've exceeded an array bound and it's specifically the second argument to the function which is Z, one would presume the most likely cause is you have a loop over the size of the array that is looking at positions Z(i),Z(i+1) and when i reaches the size of the array Z, then Z(i+1) is out of bounds.
Very common logic error to make; simply make the upper loop bound one less than the array size in the particular dimension.
Of course, you need to also think about your intersection search logic -- if it is at some boundary, can you actually reach the intersection itself by the search you are using? Work out with pencil and paper what you're actually looking at and see for the end conditions...

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Graphics Object Programming finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by