Problems with fsolve error

Hello, I'm trying to run this code:

%% BEMT Stepnieski
clear
close all
clc
%%
R=5.5;
n=100;
drdr=1/n;
drdr=[drdr:drdr:1];
dr=R/n;
r=R*drdr;
c=0.3;
rho=1.225;
theta_0=0.1396; %collective pitch [rad]
theta_twist=0.1396; %[rad]
theta_r=theta_0 - r.*theta_twist; %[rad]
b=4; 
V_tip=205;
Omega=V_tip/R; % [rad/s]
%% Application of Momentum Theory 
%Making the two thrusts equal and using fsolve
sigma=b.*c./(pi.*R);
Vc=8; 
vr0=zeros(1,100);
phir=@(x) atan((x+Vc)./Omega.*r);
alpha_r=@(x) theta_r-atan((x+Vc)./Omega.*r);
Mr=@(x) sqrt((Omega.*r).^2+(x+Vc).^2)./340;
G=@(x) 2./pi.*acos(exp(b./2.*(1-drdr)./(drdr.*(atan((x+Vc)./Omega.*r)))));
F= @(x) 0.5.*2.*pi.*alpha_r(x).*((Omega.*r).^2).*rho.*b.*c.*dr - 4.*pi...
    .*rho.*(Vc+x).*x.*r.*dr.*G(x);
[vr,fval]=fsolve(F,vr0)

but I get the following error:

Error using trustnleqn (line 28)
Objective function is returning undefined values at initial point. FSOLVE cannot continue.

Error in fsolve (line 397)
    [x,FVAL,JACOB,EXITFLAG,OUTPUT,msgData]=...

Error in BEMT_Stepnieski (line 42)
[vr,fval]=fsolve(F,vr0)

Any idea to solve it?

Answer 1

[Stephan]

Hi,

if you call your function with the initial values

F(vr0)

you get:

>> F(vr0)

ans =

   1.0e+04 *

  Columns 1 through 3

      NaN +    NaNi      NaN +    NaNi      NaN +    NaNi

  Columns 4 through 6

      NaN +    NaNi   0.0001 + 0.0000i   0.0001 + 0.0000i

  Columns 7 through 9

   0.0000 + 0.0000i  -0.0001 + 0.0000i  -0.0003 + 0.0000i

  Columns 10 through 12

The columns 1...4 return NaN, which is the problem. If you instead use

  vr0=ones(1,100);

then you get infinite imaginary parts for the first 4 values:

>> F(vr0)

ans =

   1.0e+04 *

  Columns 1 through 3

   0.0000 -    Infi   0.0000 -    Infi   0.0001 -    Infi

  Columns 4 through 6

   0.0001 -    Infi   0.0001 - 0.0765i   0.0001 - 0.0632i

  Columns 7 through 9

  -0.0000 - 0.0537i  -0.0002 - 0.0465i  -0.0004 - 0.0410i

  Columns 10 through 12

Either you correct an issue in your function, that leads to this problem or you find a valid initial array to fix this issue.

A possible workaround is:

...
drdr=[5*drdr:drdr:1];
...
vr0=zeros(1,numel(drdr));
...

which leaves the critical values out. You have to decide if this is an acceptable way.

Best regards

Stephan