Solving four coupled differential equation usig NDsolve and getting some error

Question

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];

Answer 1

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.