# Using NDSolve for 4 coupled nonlinear PDEs

I tried to solve the following coupled PDEs using NDSolve but failed.

η = 1.61;(*Refractive index*)
c = 3 10^8;(*m/s, Speed of light*)
ℏ = 1.05 10^-34;(*J/s, Planck constant in radian*)
λ = 590 10^-9;(*m,Wavelength of optical pumping*)
B = 3.87 10^21;(*J^-1s^-2m^3,Einstein B coefficient*)
Φisc = 0.625;(*ISC yield*)

pde = {D[ρ1[z, t], t] == -(η B 0.5 /c) Ibar[z, t] ρ1[z, t] + ((η B 0.5 /c) Ibar[z, t] + (1 - Φisc)/τF) ρ2[z, t],

D[ρ2[z, t], t] == (η B 0.5 /c) Ibar[z, t] ρ1[z, t] - ((η B 0.5 /c) Ibar[z,t] + 1/τF) ρ2[z, t],

D[ρ3[z, t], t] == Φisc/τF ρ2[z, t],

D[Ibar[z, t], z] + η/c D[Ibar[z, t], t] == ((-ρ1[z, t] + ρ2[z, t]) η ℏ ω B 0.5 /c) Ibar[z, t],

ρ1[z, 0] == 10^16,
ρ2[z, 0] == ρ3[z, 0] == 0,
Ibar[0, t] == 10^16};

sol = NDSolve[pde, {ρ1, ρ2, ρ3, Ibar}, {z, 0, 4 10^-3}, {t, 0, 2 10^-6},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"}}]



The warnings are about the initial and boundary conditions although I thought I've provided all known initial and boundary conditions.

• What do you try to solve? Numbers in your model are to high for numerical method. Mar 27, 2022 at 5:22
• @AlexTrounev This is a model for calculating how many and how far optical photons can penetrate a sample containing dye molecules. The numbers of optical photons and dye molecules are indeed quite large (10^16 or even more). I adopted this model from a literature of which the authors successfully got numerical results but the calculation details were not reported. Mar 27, 2022 at 8:11
• Could you give a link to the paper with the model explanation? Mar 27, 2022 at 16:58
• @AlexTrounev Sure! The link is aip.scitation.org/doi/pdf/10.1063/…. The title of the paper is "Zero-field electron spin resonance and theoretical studies of light penetration into single crystal and polycrystalline material doped with molecules photoexcitable to the triplet state via intersystem crossing". Mar 28, 2022 at 4:24
• Thank you very much. Mar 28, 2022 at 10:25

This is an extended comment, but I am hoping that it will be, at least partially, useful.

\$Version

13.0.0 for Mac OS X ARM (64-bit)


If you have a look at this equation:

   D[Ibar[z, t], z] + η/c D[Ibar[z, t], t] == ((-ρ1[z, t] + ρ2[z, t]) η ℏ ω B 0.5 /c) Ibar[z, t],


you see that you have a $$z$$- and a $$t$$-derivative for $$Ibar$$. So, if you try to include that information in the conditions, there's some progress, see below:

pde = {
D[ρ1[z, t],
t] == -(η B 0.5/c) Ibar[z, t] ρ1[z,
t] + ((η B 0.5/c) Ibar[z,
t] + (1 - Φisc)/τF) ρ2[z, t],
D[ρ2[z, t],
t] == (η B 0.5/c) Ibar[z, t] ρ1[z,
t] - ((η B 0.5/c) Ibar[z, t] + 1/τF) ρ2[z, t],
D[ρ3[z, t], t] == Φisc/τF ρ2[z, t],
D[Ibar[z, t],
z] + η/c D[Ibar[z, t],
t] == ((-ρ1[z, t] + ρ2[z,
t]) η ℏ  ω B 0.5/c) Ibar[z, t],
ρ1[z, 0] == 10^16,
ρ2[z, 0] == 0,
ρ3[z, 0] == 0,
Ibar[0, t] == 10^16,
Ibar[z, 0] == 0};


and then

NDSolveValue[pde, {ρ1[z, t], ρ2[z, t], ρ3[z, t],
Ibar[z, t]}, {z, 0, 4 10^-3}, {t, 0, 2 10^-6},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral"}}]


you get the following errors:

and the output:

Ok, so I was not inspired enough to find an appropriate value for Ibar[z, 0] -sorry- but I tried the use of ParametricNDSolve

pde = {
D[ρ1[z, t],
t] == -(η B 0.5/c) Ibar[z, t] ρ1[z,
t] + ((η B 0.5/c) Ibar[z,
t] + (1 - Φisc)/τF) ρ2[z, t],
D[ρ2[z, t],
t] == (η B 0.5/c) Ibar[z, t] ρ1[z,
t] - ((η B 0.5/c) Ibar[z, t] + 1/τF) ρ2[z, t],
D[ρ3[z, t], t] == Φisc/τF ρ2[z, t],
D[Ibar[z, t],
z] + η/c D[Ibar[z, t],
t] == ((-ρ1[z, t] + ρ2[z,
t]) η ℏ  ω B 0.5/c) Ibar[z, t],
ρ1[z, 0] == 10^16,
ρ2[z, 0] == 0,
ρ3[z, 0] == 0,
Ibar[0, t] == 10^16,
Ibar[z, 0] == xx};


and then

family =
ParametricNDSolveValue[
pde, {ρ1[z, t], ρ2[z, t], ρ3[z, t], Ibar[z, t]}, {z,
0, 4 10^-3}, {t, 0, 2 10^-6}, xx,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral"}}]


returns

• Ah, impressive!!! Thanks for your help! I am now going to play with these values and try to find out the appropriate values of rho1 and Ibar. Mar 27, 2022 at 8:14
• @Wuritianoo Glad I was able to help :-) Good luck!!!
– bmf
Mar 27, 2022 at 8:17