# NDSolve unable to handle discontinuities

I'm trying to solve a set of ordinary non-linear differential equations but running into issues. Specifically, after supplying Mathematica the below definitions and evaluating the NDSolve command in the last code sample, I get the error

NDSolve cannot handle discontinuities where the argument of DiracDelta depends on any variable besides the temporal independent variable.

DiracDeltas appear once Mathematica performs the symbolic derivatives that act on HeavisideTheta in DiscTfJ.

I believe the error message is due to DiracDelta being evaluated before all variables were replaced by their numerical values during NDSolve's operation. I tried to prevent this by placing ?NumericQ at various places behind function arguments but without any luck. Any suggestions how to get NDSolve to solve this system of ODEs would be much appreciated!

## Update

Implementing the excellent suggestions by J.M. and xzczd NDSolve now runs until about $t = -3$ at which point it throws the error

Power::infy: Infinite expression 1/Sqrt[0.] encountered.

NDSolve::nlnum: The function value {Indeterminate,0.000687141,0.00245154 +0. I,-4.25825+0. I} is not a list of numbers with dimensions {4} at {t,g2[t],Zk[t],[Lambda]k[t],[Rho]0[t],NDSolves\$4620[t],NDSolves\$4613[t],NDSolves\$4619[t],NDSolves\$4621[t]} = {-3.0213,0.,0.999229,0.467672,2.13825,-1,1,-1,-1}.

Haven't identified the cause of this yet.

## Auxiliary functions

Custom square root with branch cut along negative real $x$-axis. Second argument $y$ decides which branch to take: upper branch if $\text{Re}(y) \geq 0$, lower branch if $\text{Re}(y)<0$.

sqrt[x_?NumericQ,y_?NumericQ]=Piecewise[{{I Sqrt[-x],Re[x]<0&&Re[y]>=0},{-I Sqrt[-x],Re[x]<0&&Re[y]<0}},Sqrt[x]]

Unprotect[Power];sqrt[x_,y_]^2:=x;Protect[Power];

Derivative[1,0][sqrt][x_,y_]=1/(2sqrt[x,y]);

Derivative[0,1][sqrt][x_,y_]=0;


Partial fraction decomposition

Denominator coefficients $\alpha_k^{\pm\pm}$

\[Alpha]k[1,1][m2t_,g2t_,z_,c_]:=1/2 (1/c+m2t/z-I g2t/z)+I sqrt[1/(c z)-1/4 (1/c-m2t/z+I g2t/z)^2,-(1/c-m2t/z)];
\[Alpha]k[1,2][m2t_,g2t_,z_,c_]:=1/2 (1/c+m2t/z+I g2t/z)-I sqrt[1/(c z)-1/4 (1/c-m2t/z-I g2t/z)^2,(1/c-m2t/z)];
\[Alpha]k[2,1][m2t_,g2t_,z_,c_]:=1/2 (1/c+m2t/z-I g2t/z)-I sqrt[1/(c z)-1/4 (1/c-m2t/z+I g2t/z)^2,-(1/c-m2t/z)];
\[Alpha]k[2,2][m2t_,g2t_,z_,c_]:=1/2 (1/c+m2t/z+I g2t/z)+I sqrt[1/(c z)-1/4 (1/c-m2t/z-I g2t/z)^2,(1/c-m2t/z)];
\[Alpha]k[1,0][m2t_,g2t_,z_,c_]:=1/2 (1/c+m2t/z)+1/2 (I sqrt[1/(c z)-1/4 (1/c-m2t/z+I g2t/z)^2,-(1/c-m2t/z)]-I sqrt[1/(c z)-1/4 (1/c-m2t/z-I g2t/z)^2,(1/c-m2t/z)]);
\[Alpha]k[2,0][m2t_,g2t_,z_,c_]:=1/2 (1/c+m2t/z)-1/2 (I sqrt[1/(c z)-1/4 (1/c-m2t/z+I g2t/z)^2,-(1/c-m2t/z)]-I sqrt[1/(c z)-1/4 (1/c-m2t/z-I g2t/z)^2,(1/c-m2t/z)]);


Numerator coefficients $\beta_k^{\pm\pm}$

\[Beta]k[1,1][m2t_,g2t_,z_,c_]:=1/(2z)+I 1/(4z) (1/c-m2t/z+I g2t/z) 1 /sqrt[1/(c z)-1/4 (1/c-m2t/z+I g2t/z)^2,-(1/c-m2t/z)];
\[Beta]k[1,2][m2t_,g2t_,z_,c_]:=1/(2z)-I 1/(4z) (1/c-m2t/z-I g2t/z) 1 /sqrt[1/(c z)-1/4 (1/c-m2t/z-I g2t/z)^2,(1/c-m2t/z)];
\[Beta]k[2,1][m2t_,g2t_,z_,c_]:=1/(2z)-I 1/(4z) (1/c-m2t/z+I g2t/z) 1 /sqrt[1/(c z)-1/4 (1/c-m2t/z+I g2t/z)^2,-(1/c-m2t/z)];
\[Beta]k[2,2][m2t_,g2t_,z_,c_]:=1/(2z)+I 1/(4z) (1/c-m2t/z-I g2t/z) 1 /sqrt[1/(c z)-1/4 (1/c-m2t/z-I g2t/z)^2,(1/c-m2t/z)];
\[Beta]k[1,0][m2t_,g2t_,z_,c_]:=1/(2z)+I 1/(8z) (1/c-m2t/z+I g2t/z) 1 /sqrt[1/(c z)-1/4 (1/c-m2t/z+I g2t/z)^2,-(1/c-m2t/z)]-I 1/(8z) (1/c-m2t/z-I g2t/z) 1 /sqrt[1/(c z)-1/4 (1/c-m2t/z-I g2t/z)^2,(1/c-m2t/z)];
\[Beta]k[2,0][m2t_,g2t_,z_,c_]:=1/(2z)-I 1/(8z) (1/c-m2t/z+I g2t/z) 1 /sqrt[1/(c z)-1/4 (1/c-m2t/z+I g2t/z)^2,-(1/c-m2t/z)]+I 1/(8z) (1/c-m2t/z-I g2t/z) 1 /sqrt[1/(c z)-1/4 (1/c-m2t/z-I g2t/z)^2,(1/c-m2t/z)];


## Threshold functions

$I_j$

ker[4]=1/(2\[Pi]^2) 1/32 Sum[\[Alpha]k[i,j][m2,g2,z,c]^2 Log[\[Alpha]k[i,j][m2,g2,z,c]],{i,2},{j,2}];

tfI[0][m2_,g2_,z_,0,c_,4]=(m2^2-g2^2)/(8z^2)+Function[4#-2m2 D[#,m2]-2g2 D[#,g2]][ker[4]];

tfI[1][m2_,g2_,z_,0,c_,4]=1/(2\[Pi]^2) 1/4 m2/z^2+Function[2D[#,m2]-2m2 D[#,{m2,2}]-2g2 D[#,m2,g2]][ker[4]];

tfI[2][m2_,g2_,z_,0,c_,4]=-(1/(2\[Pi]^2)) 1/4 1/z^2+Function[2m2 D[#,{m2,3}]+2g2 D[#,{m2,2},g2]][ker[4]];


$\text{Disc}_{q_0} \, J$

discTfJ[q0_,m1_,m2_,g1_,g2_,z1_,z2_,0,c_,4]=Function[q0 D[#,q0]+2m1 D[#,m1]+2m2 D[#,m2]+2g1 D[#,g1]+2g2 D[#,g2]][
Module[{a1,a2,b1,b2},Sum[(a1=\[Alpha]k[i,0][m1,g1,z1,c];a2=\[Alpha]k[j,0][m2,g2,z2,c];b1=\[Beta]k[i,0][m1,g1,z1,c];b2=\[Beta]k[j,0][m2,g2,z2,c];
(b1 b2)/(8 \[Pi] q0) ((q0^4-2q0^2 (a1+a2)+(a1-a2)^2)/(4q0^2))^(1/2) UnitStep[Re[q0-Re[Sqrt[a1]+Sqrt[a2]]]]),{i,2},{j,2}]]];


$\partial_{q^2} \, J$

dqTfJ[q0_,m1_,m2_,g1_,g2_,z1_,z2_,0,c_,4]=1/(2\[Pi]^2) 1/2 Re[Function[-2#-q0 D[#,q0]-2m1 D[#,m1]-2m2 D[#,m2]-2g1 D[#,g1]-2g2 D[#,g2]][
Module[{a1,a2,b1,b2},1/4 Sum[a1=\[Alpha]k[i,sig1][m1,g1,z1,c];a2=\[Alpha]k[j,sig2][m2,g2,z2,c];b1=\[Beta]k[i,sig1][m1,g1,z1,c];b2=\[Beta]k[j,sig2][m2,g2,z2,c];
b1 b2 Piecewise[{{(((a1-a2)^2-(a1+a2) q0^2)/(4 q0^4 Sqrt[-(a1-a2)^2+2 (a1+a2) q0^2-q0^4]) (ArcTan[(q0^2+a1-a2)/Sqrt[-q0^4-(a1-a2)^2+2 q0^2 (a1+a2)]]+ArcTan[(q0^2-a1+a2)/Sqrt[-q0^4-(a1-a2)^2+2 q0^2 (a1+a2)]])+1/(4 q0^2)-((a1-a2) (Log[a1]-Log[a2]))/(8 q0^4)),a1-a2!=0},{1/(4q0^2)-a1/(q0^2 Sqrt[4a1 q0^2-q0^4]) ArcTan[ q0^2/Sqrt[4a1 q0^2-q0^4]],a1-a2==0}}],{i,2},{j,2},{sig1,2},{sig2,2}]]]];

dqTfJ[0, m1_, m2_, g1_, g2_, z1_, z2_, 0, c_, 4] =
1/(2 \[Pi]^2) 1/
2 Re[Function[-2 # - 2 m1 D[#, m1] - 2 m2 D[#, m2] -
2 g1 D[#, g1] - 2 g2 D[#, g2]][
Module[{a1, a2, b1, b2},
1/4 Sum[a1 = \[Alpha]k[i, sig1][m1, g1, z1, c];
a2 = \[Alpha]k[j, sig2][m2, g2, z2, c];
b1 = \[Beta]k[i, sig1][m1, g1, z1, c];
b2 = \[Beta]k[j, sig2][m2, g2, z2, c];
b1 b2 Piecewise[{{((-a1^2 + a2^2 +
2 a1 a2 (Log[a1] - Log[a2]))/(8 (a1 - a2)^3)),
a1 - a2 != 0}, {-(1/(24 a1)), a1 - a2 == 0}}], {i, 2}, {j,
2}, {sig1, 2}, {sig2, 2}]]]];


## Flow Equations

\[Rho]0Flow[\[Rho]0_,\[Lambda]k_,\[Eta]k_,g2_,Z1_,c_,d_,N_]:=-(2+\[Eta]k)\[Rho]0-(3/2 tfI[1][2\[Rho]0 \[Lambda]k,g2,Z1,0,c,d]+(N-1)/2 tfI[1][0,0,1,0,c,d])

\[Lambda]kFlow[\[Rho]0_,\[Lambda]k_,\[Eta]k_,g2_,Z1_,c_,d_,N_]:=2\[Eta]k \[Lambda]k-\[Lambda]k^2 (9/2 tfI[2][2\[Rho]0 \[Lambda]k,g2,Z1,0,c,d]+(N-1)/2 tfI[2][0,0,1,0,c,d])

g2Flow[q0_,\[Rho]0_,\[Lambda]k_,\[Eta]k_,g2_,Z1_,c_,d_,N_]:=(\[Eta]k-2)g2-2\[Rho]0 \[Lambda]k^2 (9 discTfJ[q0,2\[Rho]0 \[Lambda]k,2\[Rho]0 \[Lambda]k,g2,g2,Z1,Z1,0,c,d]-(N-1)discTfJ[q0,0,0,0,0,1,1,0,c,d])

ZkFlow[Zk_,\[Rho]0_,\[Lambda]k_,g2_,Z1_,c_,d_]:=2Zk \[Rho]0 \[Lambda]k^2 dqTfJ[0,2\[Rho]0 \[Lambda]k,0,g2,0,Z1,1,0,c,d]


## Solving the flow equations

Module[{runner = 0, counter = 0},
With[{\[Eta]k = -Zk'[t]/Zk[t], Z1 = 1, c = 1, d = 4, N = 2},
\[Rho]0\[Lambda]k\[Eta]kg2Run[1] =
NDSolve[{\[Rho]0'[
t] == \[Rho]0Flow[\[Rho]0[t], \[Lambda]k[t], \[Eta]k, g2[t],
Z1, c, d, N],
\[Lambda]k'[
t] == \[Lambda]kFlow[\[Rho]0[t], \[Lambda]k[t], \[Eta]k,
g2[t], Z1, c, d, N],
Zk'[t] ==
ZkFlow[Zk[t], \[Rho]0[t], \[Lambda]k[t], g2[t], Z1, c, d],
g2'[t] ==
g2Flow[Sqrt[
2 \[Rho]0[t] \[Lambda]k[t]/Z1], \[Rho]0[t], \[Lambda]k[
t], \[Eta]k, g2[t], Z1, c, d, N],
\[Rho]0[0] == 0.02, \[Lambda]k[0] == 0.5, Zk[0] == 1,
g2[0] == 0}, {\[Rho]0, \[Lambda]k, Zk, g2}, {t, -10, 0},
StepMonitor :>
counter++ If[Abs[t] > runner,
Print@Chop[{counter, Round[t, 1],
E^(2 t) \[Rho]0[t], \[Lambda]k[t], Zk[t], E^(2 t) g2[t]}];
runner++]]];
Plot[E^(2 t) \[Rho]0[t] /. \[Rho]0\[Lambda]k\[Eta]kg2Run[
1], {t, -10, 0}, PlotRange -> {0, All}]
Plot[\[Lambda]k[t] /. \[Rho]0\[Lambda]k\[Eta]kg2Run[1], {t, -10, 0}]
Plot[-Zk'[t]/Zk[t] /. \[Rho]0\[Lambda]k\[Eta]kg2Run[1], {t, -10, 0}]
Plot[E^(2 t) g2[t] /. \[Rho]0\[Lambda]k\[Eta]kg2Run[1], {t, -10,
0}]] // AbsoluteTiming

• Prolly should be repeated again: keep HeavisideTheta[] to symbolics, and use UnitStep[] for numerics. – J. M. will be back soon Nov 11 '17 at 1:53
• @J.M. Thanks for the tip! I'll be keeping that in mind from now on. – Casimir Nov 11 '17 at 10:16
• I mean, personally I think definition like DiscTfJ[q0_,m1_,m2_,g1_,g2_,z1_,z2_,0,c_,4] only makes the code dirty. Back to your problem, those ArcTan terms are suspicious, do you need ArcTan[y/x] or ArcTan[x, y]? – xzczd Nov 11 '17 at 12:06
• Adding Method -> StiffnessSwitching to NDSolve results in a solution (I'm not sure if this solution is reliable, StiffnessSwitching sometimes leads to wrong solution), takes about 500 seconds on my laptop. (In v9.0.1 I also need to add Re@ to the RHS of the equation system and add ?NumericQ to those *Flow function, or NDSolve won't even start calculating, but this seems to be unnecessary for your version? ) – xzczd Nov 11 '17 at 13:41
• Method -> {FixedStep, Method -> Automatic} leads to (almost) the same solution in about 343 seconds. Something similar to NDSolveValue[{0 ==(*D[u[t],t]*)D[Sqrt[u[t]^2], t] + 1, u[0] == -1}, u, {t, 0, 4}] might have happened in your code, but it's hard to go deeper given its… complexity. – xzczd Nov 11 '17 at 14:16