# Solving four coupled differential equation usig NDsolve and getting some error

Solving four coupled differential equation and one of them Dirac delta function and getting this type of error At t == 9.769456184893095, step size is effectively zero; \ singularity or stiff system suspected. Could any one help to resolve this issue.

    g1 = 1; g = 1; K1 = 0.5; \[Omega]1 = .5; \[Alpha] = \[Pi]/3;
sol1x = NDSolve[{x1'[t] ==
g1*Sqrt[10 x1[t]]
Sin[x2[t]] (-0.5*Abs[x4[t]]^2*Cos[\[Alpha]] +
0.5*Conjugate[x3[t]]*x4[t]*Sin[\[Alpha]] +
0.5*Abs[x3[t]]^2*Cos[\[Alpha]] +
0.5*Conjugate[x4[t]]*x3[t]*Sin[\[Alpha]]) -
K1 Sin[x2[t]]*Sum [DiracDelta [(t/T) - i], {i, 1, 20}],

x2'[t] == (-g1/Sqrt[
0.4*x1[t]]) (-0.5*Abs[x4[t]]^2*Cos[\[Alpha]] +
0.5*Conjugate[x3[t]]*x4[t]*Sin[\[Alpha]] +
0.5*Abs[x3[t]]^2*Cos[\[Alpha]] +
0.5*Conjugate[x4[t]]*x3[t]*Sin[\[Alpha]]) + x1[t],

x3'[t] == ((-I*(x3[
t] (0.5 \[Omega]1 +
0.5 g Sqrt[10 x1[t]] Cos[x2[t]]*Cos[\[Alpha]]) +
0.5 x4[t] g Sqrt[10 x1[t]] Cos[x2[t]]*Sin[\[Alpha]]))/
Sqrt[(Abs[x3[t]]^2 + Abs[x4[t]]^2)]),

x4'[t] == ((-I*(x4[
t] (-0.5 \[Omega]1 -
0.5 g Sqrt[10 x1[t]] Cos[x2[t]]*Cos[\[Alpha]]) +
0.5 x3[t] g Sqrt[10 x1[t]] Cos[x2[t]]*Sin[\[Alpha]]))/
Sqrt[(Abs[x3[t]]^2 + Abs[x4[t]]^2)]),
x1[0] == 1, x2[0] == \[Pi]/2, x3[0] == 0,
x4[0] == 1}, {x1, x2, x3, x4}, {t, 0, 50}, MaxSteps -> 50000,
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 5,
PrecisionGoal -> 4];

• You need to define the constant T; you haven't done so in the code you've included. For example, when I set T = 1, your code runs fine. Apr 28, 2021 at 18:56
• Nine times out of ten, the error you've described above is because Mathematica's solution has diverged. Sometimes this is a real feature of the solutions; for example, NDSolve[{x'[t] == 1/(1 - t)^2, x[0] == 1}, x[t], {t, 0, 2}] throws the same error because the solution is $x(t) = 1/(1-t)$, which diverges at $t = 1$. So it may just be that your solution really does diverge at some value of $t$. Apr 28, 2021 at 19:01
• It looks like system unstable at K1>0.4728251. Last stable solution at K1>0.4728251, T=2 Apr 28, 2021 at 20:53
• Sorry T=2 Pi is there. Apr 29, 2021 at 4:17

We extend time interval up to 41 Pi to cover all 20 impulse in Sum[DiracDelta[(t - i T)/T], {i, 1, 20}] for T=2 Pi. It looks like this system is unstable at K1>0.0662 for T=2 Pi and {t,0,41 Pi}. Last stable solution at K1=0.06614, T=2 Pi computed with this code with option Method -> "ExplicitEuler",

g1 = 1; g = 1; K1 = .06614; \[Omega]1 = .5; \[Alpha] = \[Pi]/
3; T = 2 Pi;
sol1x = NDSolve[{x1'[t] ==
g1*Sqrt[10 x1[t]] Sin[
x2[t]] (-0.5*Abs[x4[t]]^2*Cos[\[Alpha]] +
0.5*Conjugate[x3[t]]*x4[t]*Sin[\[Alpha]] +
0.5*Abs[x3[t]]^2*Cos[\[Alpha]] +
0.5*Conjugate[x4[t]]*x3[t]*Sin[\[Alpha]]) -
K1 Sin[x2[t]]*Sum[DiracDelta[(t - i T)/T], {i, 1, 20}] ,
x2'[t] == (-g1/
Sqrt[0.4*x1[t]]) (-0.5*Abs[x4[t]]^2*Cos[\[Alpha]] +
0.5*Conjugate[x3[t]]*x4[t]*Sin[\[Alpha]] +
0.5*Abs[x3[t]]^2*Cos[\[Alpha]] +
0.5*Conjugate[x4[t]]*x3[t]*Sin[\[Alpha]]) + x1[t],
x3'[t] == ((-I*(x3[
t] (0.5 \[Omega]1 +
0.5 g Sqrt[10 x1[t]] Cos[x2[t]]*Cos[\[Alpha]]) +
0.5 x4[t] g Sqrt[10 x1[t]] Cos[x2[t]]*Sin[\[Alpha]]))/
Sqrt[(Abs[x3[t]]^2 + Abs[x4[t]]^2)]),
x4'[t] == ((-I*(x4[
t] (-0.5 \[Omega]1 -
0.5 g Sqrt[10 x1[t]] Cos[x2[t]]*Cos[\[Alpha]]) +
0.5 x3[t] g Sqrt[10 x1[t]] Cos[x2[t]]*Sin[\[Alpha]]))/
Sqrt[(Abs[x3[t]]^2 + Abs[x4[t]]^2)]), x1[0] == 1,
x2[0] == \[Pi]/2, x3[0] == 0, x4[0] == 1}, {x1, x2, x3, x4}, {t,
0, 41 Pi}, Method -> "ExplicitEuler", StartingStepSize -> 0.0002, MaxSteps -> 10^6];


Visualization of stable and unstable numerical solution

LogPlot[Evaluate[Abs[{x1[t], x2[t], x3[t], x4[t]}] /. sol1x[[1]]], {t,
0.1, 41 Pi}, PlotLegends -> {x1, x2, x3, x4}, PlotRange -> All,
PlotLabel -> Row[{"K1 = ", K1}], AxesLabel -> Automatic,
PlotRange -> All]


Note, that with explicit Euler we can compute up to message General::ovfl: Overflow occurred in computation. at t=115. We can increase or decrease StartingStepSize and rationalize to get more precise solution. For the basic case with K1=0.5 we can pass error at t == 9.769456184893095 and visualize what happened with solution after that. Figure 2 shows numerical solution computed with step size 0.001 (left) and 0.0002 (right)

From the last Figure and from the message we suppose, that there is some discontinuity due to DiracDelta. To remove this discontinuity, we integrate first equation around t= i T, and finally we have

 T = 2 Pi; Integrate[
Sin[x2[t]]*Sum[DiracDelta[(t - i T)/T], {i, 1, 20}], {t, 0, 41 Pi}]

Out[]= 2 \[Pi] (Sin[x2[2 \[Pi]]] + Sin[x2[4 \[Pi]]] +
Sin[x2[6 \[Pi]]] + Sin[x2[8 \[Pi]]] + Sin[x2[10 \[Pi]]] +
Sin[x2[12 \[Pi]]] + Sin[x2[14 \[Pi]]] + Sin[x2[16 \[Pi]]] +
Sin[x2[18 \[Pi]]] + Sin[x2[20 \[Pi]]] + Sin[x2[22 \[Pi]]] +
Sin[x2[24 \[Pi]]] + Sin[x2[26 \[Pi]]] + Sin[x2[28 \[Pi]]] +
Sin[x2[30 \[Pi]]] + Sin[x2[32 \[Pi]]] + Sin[x2[34 \[Pi]]] +
Sin[x2[36 \[Pi]]] + Sin[x2[38 \[Pi]]] + Sin[x2[40 \[Pi]]])


Therefore we can use WhenEvent instead of Sum as follows

g1 = 1; g = 1; K1 = 1/2; \[Omega]1 = 1/2; \[Alpha] = \[Pi]/3; T =
2 Pi; eps = 10^-9;
sol1x = NDSolve[{x1'[t] ==
g1*Sqrt[10 x1[t]] Sin[
x2[t]] (-1/2*Abs[x4[t]]^2*Cos[\[Alpha]] +
1/2*Conjugate[x3[t]]*x4[t]*Sin[\[Alpha]] +
1/2*Abs[x3[t]]^2*Cos[\[Alpha]] +
1/2*Conjugate[x4[t]]*x3[t]*Sin[\[Alpha]]) ,
Sqrt[2/5*x1[t]] x2'[
t] == (-g1) (-1/2*Abs[x4[t]]^2*Cos[\[Alpha]] +
1/2*Conjugate[x3[t]]*x4[t]*Sin[\[Alpha]] +
1/2*Abs[x3[t]]^2*Cos[\[Alpha]] +
1/2*Conjugate[x4[t]]*x3[t]*Sin[\[Alpha]]) + x1[t],
Sqrt[(Abs[x3[t]]^2 + Abs[x4[t]]^2)] x3'[
t] == ((-I*(x3[
t] (1/2 \[Omega]1 +
1/2 g Sqrt[10 x1[t]] Cos[x2[t]]*Cos[\[Alpha]]) +
1/2 x4[t] g Sqrt[10 x1[t]] Cos[x2[t]]*Sin[\[Alpha]]))),
Sqrt[(Abs[x3[t]]^2 + Abs[x4[t]]^2)] x4'[
t] == ((-I*(x4[
t] (-1/2 \[Omega]1 -
1/2 g Sqrt[10 x1[t]] Cos[x2[t]]*Cos[\[Alpha]]) +
1/2 x3[t] g Sqrt[10 x1[t]] Cos[x2[t]]*Sin[\[Alpha]]))),
x1[eps] == 1, x2[eps] == \[Pi]/2, x3[eps] == 0, x4[eps] == 1,
WhenEvent[Mod[t, 2 Pi] == 0,
x1[t] -> x1[t] - K1 2 Pi Sin[x2[t]]] }, {x1, x2, x3, x4}, {t, eps,
50}, Method -> "ExplicitEuler", StartingStepSize -> 0.0002,
PrecisionGoal -> 5]



Note, that this code runs up to t=54.9.

• How to resolve the unstable problem. Apr 29, 2021 at 4:41
• @user68207 Do you suppose that time should be t>40 Pi`? In a case of unstable problem what do you expected to get? Apr 29, 2021 at 10:28