Finding one real solution with vpasolve

% Parameters    
Design_speed = 66; % Knots
Disp_vol = 47;     % m3
LCG = 9;           % m
B = 5;          % m
Deadrise = 7;      % degrees

nT = 0.95;         % Transmission Efficiency 
nD = 0.5;          % Overall propulsive efficiency
rho = 1027.0;      % Water density in kg/m3
omega = (1.356*10^(-6)); % Kinematic Viscosity in m2/s
g = 9.81;          % Acceleration due to gravity in m/s2

V = Design_speed*0.5148;

%% Lift coefficient CL0
Disp_force = Disp_vol * g * rho; % displacement force
CV = V/(sqrt(g*B)); % Non dimensional speed coefficient
CLbeta = Disp_force/(0.5*rho*(V^2)*(B^2));
syms CL0_t
CL0 = vpasolve(CL0_t - 0.0065*Deadrise*CL0_t^0.6 - CLbeta == 0, CL0_t);syms lambda_t
lambda = vpasolve((LCG/B)/(0.75 - (1 / (1 / ((5.236*CV^2)/(lambda_t^2)))+2.4)) == 0, lambda_t, [0,inf]);
disp(lambda)

Im trying to use vpasolve to iterate to find a real solution to lambda, but it only give a 2x1 sym where both values are 0. Ive tried rearranging the equation, and different bounds. I'm using this:DETERMINATION OF THE OPTIMUM TRIM ANGLE OF PLANING HULLS FOR MINIMUM DRAG USING SAVITSKY METHOD as the basis for it.

How could i fix this to give a 1x1 real solution?

Answer 1

[Star Strider]

When in doubt, plot the function.

Doing that here reveals that it appears to be zero at only one point, that being 0.

% Parameters    
Design_speed = 66; % Knots
Disp_vol = 47;     % m3
LCG = 9;           % m
B = 5;          % m
Deadrise = 7;      % degrees

nT = 0.95;         % Transmission Efficiency 
nD = 0.5;          % Overall propulsive efficiency
rho = 1027.0;      % Water density in kg/m3
omega = (1.356*10^(-6)); % Kinematic Viscosity in m2/s
g = 9.81;          % Acceleration due to gravity in m/s2

V = Design_speed*0.5148;

%% Lift coefficient CL0
Disp_force = Disp_vol * g * rho; % displacement force
CV = V/(sqrt(g*B)); % Non dimensional speed coefficient
CLbeta = Disp_force/(0.5*rho*(V^2)*(B^2));
syms CL0_t
CL0 = vpasolve(CL0_t - 0.0065*Deadrise*CL0_t^0.6 - CLbeta == 0, CL0_t);
syms lambda_t
lambda = vpasolve((LCG/B)/(0.75 - (1 / (1 / ((5.236*CV^2)/(lambda_t^2)))+2.4)) == 0, lambda_t, [0,inf]);
disp(lambda)

expr = (LCG/B)/(0.75 - (1 / (1 / ((5.236*CV^2)/(lambda_t^2)))+2.4));

figure
hfp = fplot(expr);
grid
axis([-1 1  -1E-2 1E-3])
X = hfp.XData;

●

X = 1×64

-1.0000 -0.9583 -0.9091 -0.8615 -0.8182 -0.7736 -0.7273 -0.6771 -0.6364 -0.5857 -0.5455 -0.5016 -0.4545 -0.4108 -0.3636 -0.3129 -0.2727 -0.2484 -0.2240 -0.2029 -0.1818 -0.1579 -0.1340 -0.1124 -0.0909 -0.0684 -0.0571 -0.0459 -0.0344 -0.0229

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Y = hfp.YData;

Y = 1×64

-0.0144 -0.0133 -0.0119 -0.0107 -0.0097 -0.0087 -0.0077 -0.0067 -0.0059 -0.0050 -0.0043 -0.0037 -0.0030 -0.0025 -0.0019 -0.0014 -0.0011 -0.0009 -0.0007 -0.0006 -0.0005 -0.0004 -0.0003 -0.0002 -0.0001 -0.0001 -0.0000 -0.0000 -0.0000 -0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Ymax,idx] = max(Y)

Ymax = 0

idx = 32

X(idx)

ans = 0

figure
fimplicit(expr, [-1 1]*1E-3, '-pb')
grid

.

Answer 2

[Torsten]

lambda = vpasolve(LCG/B*1/lambda_t==0.75-1/(5.236*CV^2/lambda_t^2+2.4) == 0, lambda_t, [0,inf]);