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!

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Hi, welcome here. Unfortunately, you code has syntax errors. Could you fix those please. $\endgroup$ Commented Jul 31, 2019 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$ Commented Jul 31, 2019 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$ Commented Jul 31, 2019 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$ Commented Jul 31, 2019 at 8:32

1 Answer 1

Reset to default
2
$\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

Improve this answer
$\endgroup$
14
  • $\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$ Commented Jul 31, 2019 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$ Commented Jul 31, 2019 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$ Commented Jul 31, 2019 at 11:09
  • $\begingroup$ @Wuritianoo What is your version? $\endgroup$ Commented Jul 31, 2019 at 12:20
  • $\begingroup$ @AlexTrounev My version is Mathematica 11.0.1.0 Mac OS. $\endgroup$ Commented Jul 31, 2019 at 12:29

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.