1
$\begingroup$

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*)
ω = 2 π (c/λ);(*rad/s,Frequency of optical pumping in radian*)
B = 3.87 10^21;(*J^-1s^-2m^3,Einstein B coefficient*)
Φisc = 0.625;(*ISC yield*)
τF = 9 10^-9;(*s,fluorescence lifetime*)

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.

Any comments are appreciated!

$\endgroup$
5
  • $\begingroup$ What do you try to solve? Numbers in your model are to high for numerical method. $\endgroup$ Mar 27, 2022 at 5:22
  • $\begingroup$ @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. $\endgroup$
    – Wuritianoo
    Mar 27, 2022 at 8:11
  • $\begingroup$ Could you give a link to the paper with the model explanation? $\endgroup$ Mar 27, 2022 at 16:58
  • $\begingroup$ @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". $\endgroup$
    – Wuritianoo
    Mar 28, 2022 at 4:24
  • $\begingroup$ Thank you very much. $\endgroup$ Mar 28, 2022 at 10:25

1 Answer 1

2
$\begingroup$

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:

error1 error2

and the output:

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

paramout

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – Wuritianoo
    Mar 27, 2022 at 8:14
  • $\begingroup$ @Wuritianoo Glad I was able to help :-) Good luck!!! $\endgroup$
    – bmf
    Mar 27, 2022 at 8:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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