1
$\begingroup$

I wish to change the term, T'[t] to 0 when S[t]<=0 when I solve ODEs. I have used WhenEvent in my NDSolve part, but it seems not work. The code is

Nrmlphoton[t_] = 4.64903*10^22*Piecewise[{{5904.29 t + 3.13641*10^8 t^2 - 7.0716*10^12 t^3 + 5.50886*10^16 t^4 - 
 1.91492*10^20 t^5 + 2.51145*10^23 t^6, 
   0 <= t < 0.0002274269422873725`}, {-2.86583*10^6 t + 1.69452*10^10 t^2 - 4.01716*10^13 t^3 + 
 4.77113*10^16 t^4 - 2.83738*10^19 t^5 + 6.75503*10^21 t^6, 
   0.0007374377891266094` < t < 0.0009481855159674386`}, {6.63274*10^6 t - 2.17507*10^10 t^2 + 2.82046*10^13 t^3 - 
 1.81037*10^16 t^4 + 5.75833*10^18 t^5 - 7.26732*10^20 t^6,
    0.0014995956366184178` < t < 0.0016925447841110407`}, {0, 
   0.0002274269422873725` <= t <= 0.0007374377891266094`}, {0, 
   0.0009481855159674386` <= t <= 0.0014995956366184178`}, {0, 
   t >= 0.0016925447841110407`}}]

OpticCoupl = 0.5;
ISCEff = 0.625;
w = 2 Pi 1.45 10^9;
Q = 3.5 10^3;
(**Tres = Q/w;**)
(**q0 =(Exp[6.63 10^-34 1.45 10^9/(1.38 10^-23 298)]-1)^-1**)
q0 = 4277;
qthermal = 4277;
B = 4 10^-8;
Kx = 2.7 10^4;
Ky = 0.6 10^4;
Kz = 1.7 10^3;
Wxz = 1.1 10^4;
Wyz = 2.2 10^4;
Wxy = 4.4 10^3;
Px = 0.76;
Py = 0.16;
Pz = 0.08;
Ns = 2.8 10^17;
sol = NDSolve[{T'[t] == Nrmlphoton ISCEff OpticCoupl, 
   S'[t] == -T'[t] + Kx X[t] + Ky Y[t] + Kz Z[t], 
   X'[t] == 
    T'[t] Px - X[t] Kx - (X[t] - Z[t]) Wxz - (X[t] - Y[t]) Wxy - 
     B (X[t] - Z[t]) q[t], 
   Y'[t] == 
    T'[t] Py - Y[t] Ky - (Y[t] - Z[t]) Wyz - (Y[t] - X[t]) Wxy, 
   Z'[t] == 
    T'[t] Pz - Z[t] Kz + (X[t] - Z[t]) Wxz + (Y[t] - Z[t]) Wyz + 
     B (X[t] - Z[t]) q[t], 
   q'[t] == -(w/Q)  (q[t] - qthermal) + B (X[t] - Z[t]) q[t],
   q[0] == q0, X[0] == Y[0] == Z[0] == T[0] == 0, S[0] == Ns, 
   WhenEvent[S[t] <= 0, T'[t] -> 0]}, {X, Y, Z, q, S}, {t, 0, 0.002}]

I got the error:

NDSolve::svnder: The variables {T'} in (T'])[t]->0 cannot be set as state variables while solving as an ODE because some are the highest order derivatives. You may be able to use these if you solve as a system of DAEs by using Method->{"EquationSimplification"->"Residual"}.

After using the Method->{"EquationSimplification"->"Residual"}, the code cannot run at all.

NDSolve::nderr: Error test failure at t == 0.`; unable to continue.

May I ask for some help? Thank you!

$\endgroup$
  • 1
    $\begingroup$ Hi, welcome here. Unfortunately, you code has syntax errors. Could you fix those please. $\endgroup$ – user21 Jul 31 at 4:43
  • $\begingroup$ @Wuritianoo Not clear what you want? If S[t]< =0, then we assume T' [t]=0. But then again S[t]>0. In this case, we assume $T'(t)\ne 0$ or leave T'[t]=0? $\endgroup$ – Alex Trounev Jul 31 at 7:50
  • $\begingroup$ @AlexTrounev Hi Alex, thanks for your reply! In my case, I would like to let T'[t]=0 when S[t]<=0 and when S[t]>0, still keep the function T'[t] == Nrmlphoton ISCEff OpticCoupl. $\endgroup$ – Wuritianoo Jul 31 at 8:12
  • $\begingroup$ @user21 Hi, thanks for notifying me the syntax errors. I have edited the piecewise function Nrmlphoton in my original post. It should now doesn't have syntax errors. $\endgroup$ – Wuritianoo Jul 31 at 8:32
1
$\begingroup$

For this problem it is necessary to use damping near S[t] = 0. Here used function (1 + Tanh[S[t]/10^8])/2.

Nrmlphoton[t_] := 
 4.64903*10^22*
  Piecewise[{{5904.29 t + 3.13641*10^8 t^2 - 7.0716*10^12 t^3 + 
      5.50886*10^16 t^4 - 1.91492*10^20 t^5 + 2.51145*10^23 t^6, 
     0 <= t < 0.0002274269422873725`}, {-2.86583*10^6 t + 
      1.69452*10^10 t^2 - 4.01716*10^13 t^3 + 4.77113*10^16 t^4 - 
      2.83738*10^19 t^5 + 6.75503*10^21 t^6, 
     0.0007374377891266094` < t < 
      0.0009481855159674386`}, {6.63274*10^6 t - 2.17507*10^10 t^2 + 
      2.82046*10^13 t^3 - 1.81037*10^16 t^4 + 5.75833*10^18 t^5 - 
      7.26732*10^20 t^6, 
     0.0014995956366184178` < t < 0.0016925447841110407`}, {0, 
     0.0002274269422873725` <= t <= 0.0007374377891266094`}, {0, 
     0.0009481855159674386` <= t <= 0.0014995956366184178`}, {0, 
     t >= 0.0016925447841110407`}}]
OpticCoupl = 0.5;
ISCEff = 0.625;
w = 2 Pi 1.45 10^9;
Q = 3.5 10^3;
(**Tres=Q/w;**)(**q0=(Exp[6.63 10^-34 1.45 10^9/(1.38 10^-23 \
298)]-1)^-1**)q0 = 4277;
qthermal = 4277;
B = 4 10^-8;
Kx = 2.7 10^4;
Ky = 0.6 10^4;
Kz = 1.7 10^3;
Wxz = 1.1 10^4;
Wyz = 2.2 10^4;
Wxy = 4.4 10^3;
Px = 0.76;
Py = 0.16;
Pz = 0.08;
Ns = 2.8 10^17; 
sol = NDSolve[{T'[
     t] == (1 + Tanh[S[t]/10^8])/2 Nrmlphoton [t] ISCEff OpticCoupl, 
   S'[t] == -T'[t] + Kx X[t] + Ky Y[t] + Kz Z[t], 
   X'[t] == 
    T'[t] Px - X[t] Kx - (X[t] - Z[t]) Wxz - (X[t] - Y[t]) Wxy - 
     B (X[t] - Z[t]) q[t], 
   Y'[t] == 
    T'[t] Py - Y[t] Ky - (Y[t] - Z[t]) Wyz - (Y[t] - X[t]) Wxy, 
   Z'[t] == 
    T'[t] Pz - Z[t] Kz + (X[t] - Z[t]) Wxz + (Y[t] - Z[t]) Wyz + 
     B (X[t] - Z[t]) q[t], 
   q'[t] == -(w/Q) (q[t] - qthermal) + B (X[t] - Z[t]) q[t], 
   q[0] == q0, X[0] == Y[0] == Z[0] == T[0] == 0, S[0] == Ns}, {T, X, 
   Y, Z, q, S}, {t, 0, 0.002}, MaxSteps -> 10^6]
{Plot[T[t] /. sol, {t, 0, .002}, AxesLabel -> {"t", "T"}], 
 Plot[S[t] /. sol, {t, 0, .002}, AxesLabel -> {"t", "S"}], 
 Plot[X[t] /. sol, {t, 0, .002}, AxesLabel -> {"t", "X"}], 
 Plot[Y[t] /. sol, {t, 0, .002}, AxesLabel -> {"t", "Y"}], 
 Plot[Z[t] /. sol, {t, 0, .002}, AxesLabel -> {"t", "Z"}], 
 Plot[q[t] /. sol, {t, 0, .002}, AxesLabel -> {"t", "g"}, 
  PlotRange -> All]}

Figure 1

$\endgroup$
  • $\begingroup$ Hi Alex, thanks a lot for your help! Your plots are indeed what I want, however, when I ran your code, I got the error "NDSolve::ndsz: At t == 0.00007564216449103049`, step size is effectively zero; singularity or stiff system suspected." Are there any specific settings I need for running the code? Thank you! May I ask is it $\endgroup$ – Wuritianoo Jul 31 at 9:47
  • $\begingroup$ I have solved the singularity issue by changing (1 + Tanh[S[t]/10^8])/2 to (1 + Tanh[S[t]/10^9])/2. Now the code seems work. $\endgroup$ – Wuritianoo Jul 31 at 10:03
  • $\begingroup$ Given the very large values of the dependent variables, does this smooth approximation actually result in T'[t]==0 when S[t]==0? Zooming in on the first graph on {t, 0, 0.0001}, I'm not sure it does. Also, check out Plot[(1 + Tanh[S[t]/10^9])/2 /. sol, {t, 0, 0.0001}, PlotRange -> {0, 1}]. $\endgroup$ – Chris K Jul 31 at 11:09
  • $\begingroup$ @Wuritianoo What is your version? $\endgroup$ – Alex Trounev Jul 31 at 12:20
  • $\begingroup$ @AlexTrounev My version is Mathematica 11.0.1.0 Mac OS. $\endgroup$ – Wuritianoo Jul 31 at 12:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.