# Stiffness or singularity in NDSolve for a system of equations with switchpoint

I am trying to solve the following set of equations:

E0 = 2. 10^11; \[Nu]0 = 0.34;
\[Mu]0 = E0/(2 (1 + \[Nu]0));
\[Lambda]0 = (E0 \[Nu]0)/((1 + \[Nu]0) (1 - 2 \[Nu]0));
\[Mu]s = 1.; Es = E0/\[Mu]0; \[Lambda]s = E0/\[Mu]0;

\[Beta] = (0.0238 10^-3)/(6.02 10^23);
\[Alpha] = ((0.0346 10^-6) (101325))/(6.02 10^23)^2;
kB = 1.38 10^-23; \[Gamma] = 1.;
pg[n_, R_, T_] := (n kB T)/((4/3) \[Pi] R^3 -
n \[Beta]) - (\[Alpha] n^2)/((4/3) \[Pi] R^3)^2;
p[n_, R_, T_] := pg[n, R, T] - (2 \[Gamma])/R;
ps[n_, R_, T_] := p[n, R, T]/\[Mu]0;

Tp = 1. 10^5; \[Tau]0 = 10^-12; \[Alpha]T = 0.03;
b = 0.387 10^-9;
sy[R_] := Sqrt[3]/2 b/R;
py[R_] := b/R;
q[s_, \[Rho]_, T_] := Sqrt[\[Rho]] Exp[-(Tp/T) Exp[-(s/(\[Alpha]T \[Mu]s Sqrt[\[Rho]]))]];
\[Lambda][\[Nu]_] := If[\[Nu] > 1, \[Nu]/(\[Nu] - 1) Log[\[Nu]] - 1, 0];
K\[Rho] = 10000;
\[Nu]bar[\[Sigma]_, p_, R_] := Max[1, CubeRoot[(\[Sigma] + p)/py[R]]];
RdotoverR1[\[Sigma]_, p_, \[Sigma]dot_, pdot_] := (4 \[Mu]s)/(4 \[Mu]s -
3 (\[Sigma] + p)) ((1 - 2 \[Nu]0)/Es \[Sigma]dot + (\[Sigma]dot +
pdot)/(4 \[Mu]s));
RdotoverR2[\[Sigma]_, p_, \[Sigma]dot_, pdot_, R_, \[Rho]_, T_] := (1 - (3 (\[Sigma] + p))/(
4 \[Mu]s (\[Nu]bar[\[Sigma], p, R])^3) + (3 b)/(
4 R) (Log[\[Nu]bar[\[Sigma], p,
R]] - \[Lambda][\[Nu]bar[\[Sigma], p,
R]]))^-1 ((\[Sigma]dot + pdot)/(
4 \[Mu]s) (1/(\[Nu]bar[\[Sigma], p, R])^3 +
3/4 \[Lambda][\[Nu]bar[\[Sigma], p, R]]) + (
Sqrt[3] q[(Sqrt[3]/2) (\[Sigma] + p), \[Rho],
T])/\[Tau]0 \[Lambda][\[Nu]bar[\[Sigma], p, R]] + (
1 - 2 \[Nu]0)/Es \[Sigma]dot);
\[Rho]dot[\[Sigma]_, p_, R_, \[Rho]_, T_] := (
2 K\[Rho])/\[Mu]s (\[Sigma] +
p) q[(Sqrt[3]/2) (\[Sigma] + p), \[Rho],
T]/\[Tau]0 (1 - \[Rho]/Exp[-4.]);

Nrate1 = 15000.; N1 = 10518.;
Ns1[t_] := N1 + Nrate1 t;
opts1 = {Method -> {"TimeIntegration" -> {"ExplicitRungeKutta",
"DifferenceOrder" -> 9}}, AccuracyGoal -> 15,
PrecisionGoal -> 15};
opts2 = {Method -> {"ExplicitRungeKutta", "StiffnessTest" -> True,
"DifferenceOrder" -> 9}, AccuracyGoal -> 15,
PrecisionGoal -> 15};
\[Epsilon] = 0.;
sol1a = NDSolve[{R'[t] ==
R[t] RdotoverR1[0, ps[Ns1[t], R[t], 300], 0,
D[ps[Ns1[t], R[t], 300], t]], R[0] == 5 10^-9}, R, {t, 0, 1},
opts1];
R1[t_] = First[R[t] /. sol1a];
switchpoint = Rationalize[t /. FindRoot[(ps[Ns1[t], R[t], 300] /. sol1a) == (py[R[t]] /. sol1a), {t, 0, 0.2}]];
sol1b = NDSolve[{R'[t] ==
R[t] RdotoverR2[0, ps[Ns1[t], R[t], 300], 0,
D[ps[Ns1[t], R[t], 300], t], R[t], \[Rho][t], 300], \[Rho]'[t] == \[Rho]dot[0, ps[Ns1[t], R[t], 300], R[t], \[Rho][t], 300],
R[switchpoint + \[Epsilon]] ==
R1[switchpoint + \[Epsilon]], \[Rho][switchpoint + \[Epsilon]] == 0.01 Exp[-4.]}, {R, \[Rho]}, {t, switchpoint + \[Epsilon], 10}, opts2];


These equations describe the expansion of a gas bubble in a solid; there is an onset of plasticity when the effective pressure reaches some critical value, hence the two-piece solution. But the solution is expected to be continuous there, if not smooth.

As it turns out, I end up with this error at the switching point:

NDSolve::ndsz: At t == 0.2191262924141748, step size is effectively zero; singularity or stiff system suspected.


I tried basically all tricks and methods that I could think of and could find in other discussions, including different StiffnessTest settings, solution methods (including Shooting), but to no avail. I did not expect a singularity at a switch point, and while the second set of equations may be stiff there, the observation that no other method has worked thus far does make me wonder. How should I proceed?

• Just rename your function called p to ps. You defined it as p but in the code you wrote ps. There is no ps function anywhere in your code. Commented Sep 26, 2023 at 20:18
• Yes, you're absolutely right. This was cut out from my research code and I forgot the one crucial line that relates p and ps. I have put that back in. Commented Sep 26, 2023 at 20:37
• You have something called \[Lambda] but it has no definition and no numerical value anywhere. Look at your RdotoverR2 function, It uses \[Lambda] as function \[Lambda][\[Nu]bar[\[Sigma], p, R]]` Commented Sep 26, 2023 at 20:49
• You're right. I missed that one too in my copy-and-paste. It's now up here. Commented Sep 26, 2023 at 20:52
• OK, with your latest code, I get 1/0 error. V 13.3.1 Screen shot !Mathematica graphics Commented Sep 26, 2023 at 20:59