2
$\begingroup$

I'm trying to solve:

$\partial _{z} U(z,t) = i \sqrt{d} P(z,t)$

$\partial _{t} P(z,t) = -P(z, t) + i \sqrt{d} U(z,t) + i \Omega_{c}(t)S(z, t)$

$\partial _{t} S(z,t) = i \Omega_{c}(t)P(z, t)$

For the initial conditions

$P(z, 0) = 0$

$S(z, 0) = 0$

$U(0, t) = A \exp\Bigg(-4 \ln(2) \Bigg(\frac{t-t_{0}}{\tau}\Bigg)\Bigg) \,,$

where A, $\tau$ and d are all constants.

As I understand it, NDSolve cannot deal with this problem due to there being a derivative of only one dimension on each of the equations. I've attempted to solve using the method of lines, as shown in answers such as this: NDSolve:Coupled PDE's, initial-boundary value problem: unreasonable "insufficient number of boundary conditions" error. However while the solver shows no errors, plotting the result just produces a blank box.

My code is shown below. I have removed the part which computes the constants to make it more legible. Capital "I" below refers to the complex number i.

d = 84.9601 tau = 0.1 OmegaC = 27.7259 pulseTime = 2 Pi/OmegaC storageTime = 0.2; twriteon = t0 - pulseTime/2; twriteoff = t0 + pulseTime/2; treadon = t0 + pulseTime/2 + storageTime; treadoff = t0 + pulseTime/2 + storageTime + pulseTime; Clear[OmegaFunc] OmegaFunc[t_?NumericQ /; t < twriteon] = 0; OmegaFunc[t_?NumericQ /; twriteoff >= t >= twriteon] = OmegaC; OmegaFunc[t_?NumericQ /; treadon > t > twriteoff] = 0; OmegaFunc[t_?NumericQ /; treadoff >= t >= treadon] = OmegaC; OmegaFunc[t_?NumericQ /; t > treadoff] = 0; (*Plot[OmegaFunc[t], {t, 0, tend}, PlotRange -> All]*) (*Initial conditions*) p0 [z_] = 0; s0[z_] = 0; A = 9.39437 u0[t_] = A* Exp[-4*Log[2]*((t - t0 )/tau)^2]; (*Plot[u0[t], {t, 0, tend}, PlotRange -> All]*) (*Define arrays in z to discretize problem in z*) zmax = 1; n = 4; h = zmax/n; P[t_] = Table[p[i][t], {i, 1, n}]; S[t_] = Table[s[i][t], {i, 1, n}]; integrateP = Join[{p[1][t]}, Table[p[i - 2][t] + p[i - 1][t], {i, 3, n}]]; U[t_] = Join[{u0[t]}, integrateP]; (*Construct equations*) eqP = Thread[ D[P[t], t] == -P[t] + I Sqrt[d] U[t] + I OmegaFunc[t] S[t]]; eqS = Thread[D[S[t], t] == I OmegaFunc[t] P[t]]; initP = Thread[P[0] == Table[p0[(i - 1) h], {i, 1, n}]]; initS = Thread[S[0] == Table[s0[(i - 1) h], {i, 1, n}]]; (*Solve*) lines = NDSolve[{eqP, eqS, initP, initS}, {P[t], S[t]}, {t, 0, tend}]; (*Plot*) ztab = Table[(i - 1) h, {i, 1, n}]; ParametricPlot3D[ Evaluate@Thread[{ztab, t, lines[[1, 1]]}], {t, 0, tend}, PlotRange -> All, AxesLabel -> {"z", "t", "P"}, BoxRatios -> {2, 2, 1}, ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}] ParametricPlot3D[ Evaluate@Thread[{ztab, t, lines[[1, 2]]}], {t, 0, tend}, PlotRange -> All, AxesLabel -> {"z", "t", "S"}, BoxRatios -> {2, 2, 1}, ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

I've shown the result of plotting P(z, t) below. S(z, t) is the same.

1

Could anyone help with what I'm doing wrong? Thanks in advance.

Improve this question
$\endgroup$
1
  • $\begingroup$ z is over the range [0, 1] and t is [0, 1.4] $\endgroup$ Commented Aug 6, 2020 at 14:06

1 Answer 1

Reset to default
2
$\begingroup$

We can use first order approximation on z as follows

d = 84.9601; tau = 0.1; t0 = .1; OmegaC = 27.7259; pulseTime = 2 Pi/OmegaC; storageTime = 0.2; twriteon = t0 - pulseTime/2; twriteoff = t0 + pulseTime/2; treadon = t0 + pulseTime/2 + storageTime; treadoff = t0 + pulseTime/2 + storageTime + pulseTime; Clear[OmegaFunc] OmegaFunc[t_] := Piecewise[{{0, t < twriteon}, {OmegaC, twriteoff >= t >= twriteon}, {0, treadon > t > twriteoff}, {OmegaC, treadoff >= t >= treadon}, {0, t > treadoff}}] (*Plot[OmegaFunc[t],{t,0,tend},PlotRange\[Rule]All]*) (*Initial conditions*) p0[z_] = 0; s0[z_] = 0; A = 9.39437; u0[t_] = A*Exp[-4*Log[2]*((t - t0)/tau)^2]; (*Plot[u0[t],{t,0,tend},PlotRange\[Rule]All]*) (*Define arrays in z to discretize problem in z*) zmax = 1; n = 41; h = zmax/n; P[t_] = Table[p[i][t], {i, 1, n}]; S[t_] = Table[s[i][t], {i, 1, n}]; (*Construct equations*) eqP = Table[ D[p[i][t], t] == -p[i][t] + I Sqrt[d] (u0[t] + h Sum[p[j][t], {j, 1, i}]) + I OmegaFunc[t] s[i][t], {i, n}]; eqS = Table[D[s[i][t], t] == I OmegaFunc[t] p[i][t], {i, n}]; initP = Table[p[i][0] == 0, {i, 1, n}]; initS = Table[s[i][0] == 0, {i, 1, n}]; (*Solve*) tend = 1.4; var = Join[Table[p[i], {i, n}], Table[s[i], {i, n}]]; sols = NDSolve[{eqP, eqS, initP, initS}, var, {t, 0, tend}]; (*Plot*)

Since functions P, S, U are complex we can visualize Re,Im, Abs as

ztab = Table[(i - 1) h, {i, 1, n}]; U = Join[{u0[t]}, h Table[Sum[p[j][t], {j, i}], {i, n}]]; lst1 = Flatten[ Table[{ztab[[i]], t, Re[p[i][t]] /. First[sols]}, {i, n}, {t, 0, tend, .01 tend}], 1]; lst2 = Flatten[ Table[{ztab[[i]], t, Re[s[i][t]] /. First[sols]}, {i, n}, {t, 0, tend}], 1]; lst3 = Flatten[Table[{ztab[[i]], t, Re[U[[i]]] /. First[sols]}, {i, n}, {t, 0, tend}], 1]; {ListPlot3D[lst1, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "P"}, PlotRange -> All, Mesh -> None], ListPlot3D[lst2, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "S"}, PlotRange -> All, Mesh -> None], ListPlot3D[lst3, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "U"}, PlotRange -> All, Mesh -> None]} lst11 = Flatten[ Table[{ztab[[i]], t, Im[p[i][t]] /. First[sols]}, {i, n}, {t, 0, tend, .01 tend}], 1]; lst21 = Flatten[ Table[{ztab[[i]], t, Im[s[i][t]] /. First[sols]}, {i, n}, {t, 0, tend}], 1]; lst31 = Flatten[Table[{ztab[[i]], t, Im[U[[i]]] /. First[sols]}, {i, n}, {t, 0, tend}], 1]; {ListPlot3D[lst11, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "P"}, PlotRange -> All, Mesh -> None], ListPlot3D[lst21, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "S"}, PlotRange -> All, Mesh -> None], ListPlot3D[lst31, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "U"}, PlotRange -> All, Mesh -> None]} lst12 = Flatten[ Table[{ztab[[i]], t, Abs[p[i][t]] /. First[sols]}, {i, n}, {t, 0, tend, .01 tend}], 1]; lst22 = Flatten[ Table[{ztab[[i]], t, Abs[s[i][t]] /. First[sols]}, {i, n}, {t, 0, tend}], 1]; lst32 = Flatten[Table[{ztab[[i]], t, Im[U[[i]]] /. First[sols]}, {i, n}, {t, 0, tend}], 1]; {ListPlot3D[lst12, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "|P|"}, PlotRange -> All, Mesh -> None], ListPlot3D[lst22, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "|S|"}, PlotRange -> All, Mesh -> None], ListPlot3D[lst32, ColorFunction -> "Rainbow", AxesLabel -> {"z", "t", "|U|"}, PlotRange -> All, Mesh -> None]}

Figure 1

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.