Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. when switching to stiff solvers

Question

Yi-xiao Liu 2021-3-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/780627-warning-matrix-is-singular-close-to-singular-or-badly-scaled-results-may-be-inaccurate-rcond-n

回答： Yi-xiao Liu 2021-3-30

I am solving a system of ODE to model electrospinning, which could be described as n Viscoelastically connected point charge

    y0=[r(LeadingBeadidx:LastBeadidx+1,:),v(LeadingBeadidx:LastBeadidx+1,:),...
        sigmau(LeadingBeadidx:LastBeadidx+1)];
    y0=reshape(y0,[],1);
    %ODE45 expects a column vector
    [t,YSol]=ode45(@(t,y) MasterODE(t,y,m(LeadingBeadidx:LastBeadidx+1),...
        e(LeadingBeadidx:LastBeadidx+1),alpha,E,G,mu,InitSegV),...
        [0,dt],y0,ODEoptions);
    YSol=reshape(YSol(t==dt,:),[],7);
function dVardt= MasterODE(t,Var,m,e,alpha,E,G,mu,InitSegV) %#ok<INUSL>
%Mater ODE that combines motion(newton law) and sigma(maxwell materials for
%all beads
GPUCompute=numel(Var)>=700 && gpuDeviceCount>=1;
if GPUCompute
    try
        %transfer larger data to GPU, if present
        Var=gpuArray(Var);
        m=gpuArray(m);
        e=gpuArray(e);
    catch
    end
end
Var=reshape(Var,[],7);
r=Var(:,1:3);
v=Var(:,4:6);
sigmau=Var(:,7);
%========================================================================
%calculate required parameters from r and v for force calculation
%========================================================================
f=Viscoelastic(au,sigmau,dispu,lu)-Viscoelastic(ad,sigmad,dispd,ld);
%Calculate f, starting with viscoelastic forces
f=f+SurfaceTension(au,ad,r,dispcross,dispu,dispd,alpha);
%adding surface tension
f=f+ExternalField(e,E);
%Adding external electric field
f=f+InterElectric(e,r);
%Adding inter beads electrostatic forces
%========================================================================
drdt=v;
%Derivative of position is velocity
dvdt=f./m;
%f=ma,a=dvdt
dsigmaudt=(G./lu).*dludt-(G./mu).*sigmau;
%Derivatives of sigma
dVardt=[drdt,dvdt,dsigmaudt];
dVardt(isnan(dVardt))=0;
%Assemble derivatives
dVardt=reshape(dVardt,[],1);%ODE45 expects column vector
if GPUCompute
    dVardt=gather(dVardt);% Return array to workspace
end
end

where r and v are n-by-3 vectors for the position and velocity of n points, sigma is a n-by-1 vector of stress between points.

This system could be solved by ode45 reasonably fast when the n is small. However when n is large (like >1000) it takes forever, even for a small time span. Feeling that the system may become stiff, I was tempted to try some stiff solvers like ode15s, but got greeted by the following error:

...
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
Warning: Failure at t=0.000000e+00.  Unable to meet integration tolerances without reducing the step size below the smallest
value allowed (7.905050e-323) at time t. 