# Solving this system of differential equation with NDSolve

I'm trying to solve this system of differential equation (the ones in yellow markers are just constants): The initial conditions for this system are n(0) = NN, na(0) = Na, n'(0) = q(0) = 0, q'(0) = some_constants, na'(0) = Ip(0)×Some_other_constants (shown in my code below).

FWHM = 9;
sd = 2.355*FWHM;
\[Sigma]p = 1.15*10^-21; (* m^2 *)
\[Sigma]pa = 3.8*10^-21; (* m^2 *)
\[Sigma]e11 = 1.1*10^-20;  (* m^2 *)
\[Sigma]111 = 0.4*10^-20;  (* m^2 *)
\[Sigma]aa1 = 4.1*10^-20;  (* m^2 *)
\[Sigma]ea1 = 0.9*10^-20;  (* m^2 *)
\[Sigma]1a1 = 1.3*10^-20;  (* m^2 *)
\[Sigma]eaa = 3.05*10^-20;  (* m^2 *)
\[Sigma]1aa = 2.0*10^-20;  (* m^2 *)

NN = 9*10^27; (* m^-3 *)
Na = 1.08*10^26; (* m^-3 *)

c = 3*10^8; (* m/s*)
\[Eta] = 1.32;
\[CapitalOmega] = 5.2*10^-8;
\[CapitalOmega]a = 2.5*10^-6;
L = 0.003; (*m*)
V = 0.38;
\[Tau] = 4*10^-9 ;(*s*)
\[Tau]a = 4*10^-9; (*s*)

Ip [t] ==  2.5 E^(-(t - 8)^2/(2 sd))*10^28;

sol = NDSolve[{
n'[t] == Ip[t] \[Sigma]p (NN - n[t]) - (\[Sigma]e11 c)/\[Eta] n[t] q[t] -
n[t]/\[Tau],
q'[t] == ((\[Sigma]e11 - \[Sigma]111) c)/\[Eta] n[t] q[t] -
q[t]/\[Tau]c[t] + (\[CapitalOmega] n[t])/\[Tau] -
(\[Sigma]aa1 - c)/\[Eta] na[t]q[t] + (\[Sigma]ea1*c)/\[Eta] (Na - na[t]) q[t] -
(\[Sigma]aa1 c)/\[Eta] (Na - na[t]) q[t],
na'[t] == -Ip[t] \[Sigma]pa na[t] + (\[Sigma]ea1 c)/\[Eta] (Na - na[t]) q[t] +
(Na - na[t])/\[Tau]a - (\[Sigma]aa1 c)/\[Eta] na[t] q[t] +
(2 \[CapitalOmega]a)/(\[Tau]a L (\[Sigma]eaa - \[Sigma]1aa))*
((exp ((\[Sigma]eaa - \[Sigma]1aa) (Na - na[t]) L)) - 1)  ,
\[Tau]c [t_] == (\[Eta] L^3)/(8 c \[Pi]^2) (n[t] (\[Sigma]e11 - \[Sigma]111) V)^2,
Ip [t_] ==  2.5 E^(-(t - 8)^2/(2 sd))*10^28,

n == NN , q == 0, na == Na,
n' == 0, q' == (\[CapitalOmega] NN)/\[Tau],
na' == -Ip \[Sigma]pa Na,
\[Tau]c == (\[Eta] L^3)/(8 c \[Pi]^2) (NN (\[Sigma]e11 - \[Sigma]111) V)^2},
{n, q, na}, {t, 0, 17}];

Plot[Evaluate[{n[t], na[t], q[t], \[Tau][t]} /. s], {t, 0, 17},
PlotPoints -> (t, 0, 1000)]


I've tried to write my code as shown above, however I got an error message

NDSolve::idelay: Initial history needs to be specified for all variables for delay-differential equations.

I'm not familiar with this topic and I personally think that I am not intending to solve a 'delay' equation. Can anyone pointing out my mistakes? Thank you.

• Please check your equation. What is exp()? It should be Exp[]. And should be Ip [t_] = 2.5 E^(-(t - 8)^2/(2 sd))*10^28 Jun 12, 2022 at 23:35
• \[Tau]c [t_] needs to come outside of NDSolve and have correct function syntax. Same with Ip [t]. But even with all these fixes I'm still getting The number of constraints 6 is not equal to the total differential order of the system plus the number of discrete variables 3 and that means there are still more errors in this.
– Bill
Jun 12, 2022 at 23:50

I fixed the exp that @cvgmt noted, pulled τc and Ip out of NDSolve, and got rid of the extra initial conditions.

τc[t] := (η L^3)/(8 c π^2) (n[t] (σe11 - σ111) V)^2;
Ip[t] := 2.5 E^(-(t - 8)^2/(2 sd))*10^28;

eqns = {n'[t] == Ip[t] σp (NN - n[t]) - (σe11 c)/η n[t] q[t] - n[t]/τ,
q'[t] == ((σe11 - σ111) c)/η n[t] q[t] - q[t]/τc[t] + (Ω n[t])/τ
- (σaa1 - c)/η na[t] q[t] + (σea1*c)/η (Na - na[t]) q[t]
- (σaa1 c)/η (Na - na[t]) q[t],
na'[t] == -Ip[t] σpa na[t] + (σea1 c)/η (Na - na[t]) q[t] + (Na - na[t])/τa
- (σaa1 c)/η na[t] q[t] + (2 Ωa)/(τa L (σeaa - σ1aa))*(Exp[(σeaa - σ1aa) (Na - na[t]) L] - 1)};

sol = NDSolve[{eqns, n == NN, q == 0, na == Na}, {n, q, na}, {t, 0, 17}];


which results in

NDSolve::nderr -- Error test failure at t == 1.740427851496955*^-8; unable to continue.

To get a hint why that might happen, I evaluated the time-derivatives at t==0:

eqns /. t -> 0 /. {n -> NN, q -> 0, na -> Na}
(* {n' == -2.25*10^36, q' == 1.17*10^29, na' == -2.26701*10^33} *)


Those are some absolutely huge numbers, so it's not surprising that NDSolve can't deal with this over a time range of t==0 to t==17. Are you sure those parameter values are correct? (I guess if the variables are also of large magnitude, it might OK?)

ps. out of curiosity, what do the equations model?

• Also, I just noticed that the fourth term in q'[t] should have \[Sigma]aa1*c not \[Sigma]aa1- c. Fixing that gets NDSolve to t==1.9939`. So go over your whole problem to look for other typos if the parameters are in fact correct. Jun 13, 2022 at 2:44