No solution found. fsolve stopped because the problem appears to be locally singular.
3 views (last 30 days)
Show older comments
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
4 Comments
Torsten
on 29 Nov 2023
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.
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!