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

Question

Athanasios 2023-11-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2053847-no-solution-found-fsolve-stopped-because-the-problem-appears-to-be-locally-singular

评论： Athanasios 2023-11-30

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 个评论
显示 2更早的评论隐藏 2更早的评论

Athanasios 2023-11-29

编辑：Walter Roberson 2023-11-29

在 MATLAB Online 中打开

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