Matrix Near Singularity Warning During the ODEs Solving Procedure

I am working with multi rigid body systems simulation utilizing Runge-Kutta (RK) numerical integration method. In the every RK iteration step, there is a linear equations to solve in order to obtain the state's derivative. The ODEs is illustrated as follows：

function Ydot = ODEs(t,Y,params)
    % A simple demonstration of ODEs
    % Y=[q;qdot] where q: generalized positions
    %               qdot: generalized velocities
    numq = length(params.q0); % number of generalized coordinates
    q = Y(1:numq); % 
    qdot = Y(numq+1:end);
    
    Ydot = zeros(2*numq,1); % preallocation Ydot
    Ydot(1:numq) = qdot;
    LHS = A(q,qdot); % get linear equations' left hand side matrix
    RHS = b(q,qdot); % get linear equations' right hand side vector (Ax=b)
    Ydot(numq+1:end) = LHS\RHS; % obtain generalized acceleration which produce warning!
end

function rhs = b(q,qdot)
    % function for obtaining right hand side vector
    % A simple illustration
    % 
    delta = getDelta(q,qdot); % obtain delta value. About 1e-6 orders of magnitude
    rhs = max(0,delta).^1.5*other; % calculate rhs 
end

However, when I simulate the ODEs with the exponent 1.5 setted as shown in max(0,delta).^1.5. The matrix near singularity warning occurs in the solving procedure of LHS\RHS. When I changed the exponent 1.5 to 2, the simulation process is normal and there is no warning.

I suspect the warning is caused by the smaller exponent of max(0,delta), since the value of delta is very small and the larger exponent brings the max(0,delta) closed to 0. The smaller value of max(0,delta).^2 mitigate the problem of badly number scaling.

I want to simulate under different exponent (e.g., 1~3), can anyone give me some advice how to fix this warning when the exponent is small

Answer 1

[Chuguang Pan]

Through stepwise debugging, I find that the badly scaling issue was induced by the stiffness value K, which is a large value. When I decrease the value of K, the ODEs are solved normally without warning. I think the large value of stiffness results in the ODEs to be stiff.