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

4 visualizaciones (últimos 30 días)
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.
So my questions are:
1.What is the cause of this error? How do I solve it?
2.I could supply part of the jacobian, or at least a jPattern. Will this help?
3. Does this error means that this is not a stiff problem and not a good fit for stiff solvers? How does one determine when the system is stiff?
Thanks!

Respuesta aceptada

Yi-xiao Liu
Yi-xiao Liu el 30 de Mzo. de 2021
OK, I am going to answer my own question here.
After providing Jpattern, the ode15s works normally, but does not display any performance advantages over ode45. But at least Jpattern makes it work.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Productos


Versión

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by