No solution found. fsolve stopped because the problem appears to be locally singular.

I'm using this code in order to model a bearing with 5 degrees of freedom. I'm having the 5 nonlinear differential equations that describe the model and I'm using the Newmark method, in order to transform them into nonlinear algebraic equations. At this phase I'm using the 'fsolve' solver so that I can solve the nonlinear algebraic system for each time step. Unfortunately, I get the message 'No solution found. fsolve stopped because the problem appears to be locally singular.' I have seen other answers in the forum with the same topic, but still I have no clue what's wrong. I have also used 'vpasolve' function but I'm getting empty double column vector. Thanks to everyone in advanced.

% Define the bearing parameters
Dp = 98.5/1000; % m
Db = 15/1000; % m
Nball = 10; %
Di = 83.5/1000; % m
Do = 113.5/1000; % m
Cr = 0.05/1000; % m
a = 0; % degrees/rads
kb = 1.89e8; % N/m
kp = 1.51e7; % N/m
ks = 5.24e4; % N/m
kr = 1.88e9; % N/m
cp = 2310.68; % Ns/m
cs = 3376.84; % Ns/m
cr = 10424.46; % Ns/m
mp = 4.34; % kg
ms = 2.12; % kg
mr = 1.5; % kg
ws_rpm = 1000; % shaft speed in rpm
ws = 2*pi*ws_rpm/60; % convert shaft speed from rpm to rad/s 
wc = (1 - Db*cos(a)/Dp)*ws/2; % cage speed in rad/s
theta0 = 0; % initial angular position of the first rolling element relative to x coordinates in rad
M = [ms 0 0 0 0
    0 ms 0 0 0
    0 0 mp 0 0 
    0 0 0 mp 0
    0 0 0 0 mr];
C = [cs   0    0      0     0
    0   cs    0      0     0
    0    0   cp      0     0 
    0    0    0  cp+cr   -cr
    0    0    0    -cr    cr];
K = [ks   0    0      0     0
     0   ks    0      0     0
     0    0   kp      0     0 
     0    0    0  kp+kr   -kr
     0    0    0    -kr    kr];
tf = 2; % 2 seconds
a_newmark = 0.25;
b_newmark = 0.5;
n = size(M,1); 
h = 1e-3;
E1 = (1/h^2)*M + (b_newmark/h)*C + a_newmark*K;
E2 = (1/h^2)*M + (b_newmark/h)*C;
E3 = (1/h)*M + (b_newmark-a_newmark)*C;
E4 = (0.5 - a_newmark)*M + 0.5*h*(b_newmark - 2*a_newmark)*C;
t = 0:h:tf;
x = zeros(n,length(t));
dx = zeros(n,length(t));
ddx = zeros(n,length(t));
Fr = 3000;
F = [0; Fr; 0; 0; 0]; % external Force vector 
%%%%%%%%% NEWMARK %%%%%%%%%
for i = 2:length(t)
    disp(t(i))
    % Define angular positions of the balls
    phi = zeros(Nball,1);
    for j = 1:Nball
        phi(j,1) = mod(theta0 + wc*t(i) + 2*pi()*(j-1)/Nball, 2*pi()); % angular positions of balls
    end
    
    % Define nonlinear forces fx and fy
    syms xs ys xp yp yb
    delta = (xs - xp)*sin(phi-3*pi()/2) + (ys - yp)*cos(phi-3*pi()/2) - cr/2;
    n= 10/9; 
    delta_stin_n_cos= delta.^n.*cos(phi);
    delta_stin_n_sin= delta.^n.*sin(phi);
    fx = kb * sum(delta_stin_n_cos);
    fy = kb * sum(delta_stin_n_sin);
    % Convert symbolic expressions to double-precision functions
    fx_numeric = matlabFunction(fx, 'Vars', {xs, ys, xp, yp, yb});
    fy_numeric = matlabFunction(fy, 'Vars', {xs, ys, xp, yp, yb});
    % Solve nonlinear algebraic equation 
    vars = [xs; ys; xp; yp; yb];
    eqn_function = @(vars) [E1(1, :) * vars + a_newmark * fx_numeric(vars(1), vars(2), vars(3), vars(4), vars(5)) -     E2(1, :) * x(:, i-1) - E3(1, :) * dx(:, i-1) - E4(1, :) * ddx(:, i-1) - a_newmark * F(1, 1);
        E1(2, :) * vars + a_newmark * fy_numeric(vars(1), vars(2), vars(3), vars(4), vars(5)) - E2(2, :) * x(:, i-1) - E3(2, :) * dx(:, i-1) - E4(2, :) * ddx(:, i-1) - a_newmark * F(2, 1);
        E1(3, :) * vars - a_newmark * fx_numeric(vars(1), vars(2), vars(3), vars(4), vars(5)) - E2(3, :) * x(:, i-1) - E3(3, :) * dx(:, i-1) - E4(3, :) * ddx(:, i-1) - a_newmark * F(3, 1);
        E1(4, :) * vars - a_newmark * fy_numeric(vars(1), vars(2), vars(3), vars(4), vars(5)) - E2(4, :) * x(:, i-1) - E3(4, :) * dx(:, i-1) - E4(4, :) * ddx(:, i-1) - a_newmark * F(4, 1);
        E1(5, :) * vars + a_newmark * 0 - E2(5, :) * x(:, i-1) - E3(5, :) * dx(:, i-1) - E4(5, :) * ddx(:, i-1) - a_newmark * F(5, 1)];
    % Define your options
    options = optimoptions('fsolve', 'FunctionTolerance', 1e-6, 'StepTolerance', 1e-8);
    % Initial guess for the solution
    initial_guess = [0;0;0;0;0];
    % Solve the system numerically using fsolve with the specified options
    sol = fsolve(eqn_function, initial_guess, options);
    % Extract numerical values
    xs_value = sol(1);
    ys_value = sol(2);
    xp_value = sol(3);
    yp_value = sol(4);
    yb_value = sol(5);
    
    % ddx(:,i) = (1/(a_newmark*h^2))*(x(:,i) - x(:,i-1)) - (1/(a_newmark*h))*dx(:,i-1) + ...
    %     (1 - 1/(2*a_newmark))*ddx(:,i-1);
    % dx(:,i) = (b_newmark/(a_newmark*h))*(x(:,i) - x(:,i-1)) + (1-(b_newmark/a_newmark))*dx(:,i-1) + ...
    %     h*(1-(b_newmark/(2*a_newmark)))*ddx(:,i-1);
end

[Walter Roberson]

I turned in into an optimization of the residue of the sum of squares of the individual equations, to get an idea of how close you could get: if each individual component as perfectly solved, the contribution from that component would be 0.

I used ga with population size 2500. None of the runs came anywhere close to solving the equations; the minimum norm I saw was over 1 1/2 million.

Is it possible that your solutions are complex-valued ?

[Athanasios]

You are right beacuse I tried solving the same system with ''ode45'' function and I got complex valued solutions. I wonder how is that possible if the nonlinear differential equations represent a physical system.

Here is the algorithm if needed.

% Define the bearing parameters
Dp = 98.5/1000; % m
Db = 15/1000; % m
Nball = 2;
Di = 83.5/1000; % m
Do = 113.5/1000; % m
Cr = 0.05/1000; % m
a = 0; % degrees/rads
kb = 1.89e8; % N/m
kp = 1.51e7; % N/m
ks = 5.24e4; % N/m
kr = 1.88e9; % N/m
cp = 2310.68; % Ns/m
cs = 3376.84; % Ns/m
cr = 10424.46; % Ns/m
mp = 4.34; % kg
ms = 2.12; % kg
mr = 1.5; % kg
ws_rpm = 500; % shaft speed in rpm
ws = 2*pi*ws_rpm/60; % convert shaft speed from rpm to rad/s 
wc = (1 - Db*cos(a)/Dp)*ws/2; % cage speed in rad/s
theta0 = 0; % initial angular position of the first rolling element relative to x coordinates in rad
% Define the system of nonlinear differential equations
syms xs(t) ys(t) xp(t) yp(t) yb(t) T Y
phi = sym(zeros(Nball,1));
for j = 1:Nball
    phi(j,1) = theta0 + wc*T + 2*pi()*(j-1)/Nball;
end
delta = (xs - xp)*sin(phi-3*pi()/2) + (ys - yp)*cos(phi-3*pi()/2) - cr/2;
n= 10/9;
delta_stin_n_cos= delta.^n.*cos(phi);
delta_stin_n_sin= delta.^n.*sin(phi);
fx = kb * sum(delta_stin_n_cos);
fy = kb * sum(delta_stin_n_sin);
f1 = ms*diff(diff(xs)) + cs*diff(xs) + ks*xs + fx == 0;
f2 = ms*diff(diff(ys)) + cs*diff(ys) + ks*ys + fy == 0;
f3 = mp*diff(diff(xp)) + cp*diff(xp) + kp*xp - fx == 0;
f4 = mp*diff(diff(yp)) + (cp + cr)*diff(yp) + (kp + kr)*yp - kr*yb - cr*diff(yb) - fy == 0;
f5 = mr*diff(diff(yb)) + cr*(diff(yb) - diff(yp)) + kr*(yb - yp) == 0;
Eqns = [f1; f2; f3; f4; f5];
[DEsys,Subs] = odeToVectorField(Eqns);
DEfcn = matlabFunction(DEsys,'Vars',{T, Y});
% Set up initial conditions
initial_conditions = [0 0.02 0 0 0 0 0 0 0 0];
tspan = [0, 3];
[T,Y] = ode45(DEfcn, tspan, initial_conditions);
%% Obtain forces fx and fy
% Extract the results for each variable
xs_result = Y(:, 1);
ys_result = Y(:, 2);
xp_result = Y(:, 3);
yp_result = Y(:, 4);
% Initialize arrays to store fx and fy values
fx_values = zeros(size(T));
fy_values = zeros(size(T));
% Calculate fx and fy for each time step
for i = 1:length(T)
    t = T(i);
    xs_t = xs_result(i);
    ys_t = ys_result(i);
    xp_t = xp_result(i);
    yp_t = yp_result(i);
    
    % Update phi and delta using the current values
    phi_t = theta0 + wc*t + 2*pi()*(0:Nball-1)/Nball;
    delta_t = (xs_t - xp_t)*sin(phi_t-3*pi()/2) + (ys_t - yp_t)*cos(phi_t-3*pi()/2) - cr/2;
    delta_stin_n_cos_t = delta_t.^n.*cos(phi_t);
    delta_stin_n_sin_t = delta_t.^n.*sin(phi_t);
    
    % Calculate fx and fy for each time step
    fx_values(i) = kb * sum(delta_stin_n_cos_t);
    fy_values(i) = kb * sum(delta_stin_n_sin_t);
end

[Torsten]

I wonder how is that possible if the nonlinear differential equations represent a physical system.

Most probably, delta^n (delta_t^n) become complex-valued if delta (delta_t) become negative.

[Athanasios]

This is it! Thank you very much. I appreciate your help.