# Solving a system of delayed partial differential equations

Dear (more advanced) users of Mathematica,

I'm still a beginner and currently trying to solve a system of delayed partial differential equations. Two of them describe the evolution of of two optical fields ($s_\text{P}$ and $s_\text{St}$) along a waveguide (position z):

(1) $\frac{\partial s_\text{P}\left(z,t\right)}{\partial z}=i \kappa_\text{S} R\left(z,t\right) s_\text{St}(z,t)$, and (2) $\frac{\partial s_\text{St}\left(z,t\right)}{\partial z}=i \kappa_\text{S} R^*\left(z,t\right) s_\text{P}\left(z,t\right)$.

These equation are coupled by the amplitude $R\left(z,t\right)$ of vibrations driven by the optical fields, which itself depends on time:

(3) $\frac{\partial R\left(z,t\right)}{\partial t} + \frac{\Gamma}{2} R\left(z,t\right) = i \kappa_\text{R} s_\text{P}\left(z,t\right) s_\text{St}^* \left(z,t\right)$.

I can obtain the steady-state solution ($\partial R/\partial t =0$) very easily by putting (3) into (1),(2) and applying NDSolve and the boundary conditions ($s_0\left(0,t\right) = \sqrt{0.9}$ and $s_\text{-1}(0,t)=\sqrt{0.1}$) to these two equations. But what I'd be interested in is the build-up of R over time. If I simply integrate (3) with $R(z,0) = 0$, I end up with an "delayed PDE" error.

Does anyone have an idea how I could proceed?

EDIT: This is my kind of naive code to solve this problem:

    (*DEFINE ALL THE CONSTANTS*)

Clear[ΩR, αS, κR, κS, γR, sar, c0];
ω0 = 2*Pi*2.99792458*10^8/(1.55*10^(-6));   (*optical frequency in Hz*)
c0 = 2.99792458*10^8;                       (*speed of light in vacuum*)
ΩR = 2*Pi*5.611*10^6;                       (*mechanical frequency in Hz*)
γR = 2*Pi*400;                              (*damping rate/linewidth in Hz*)
κS = -219.666*0.1;                          (*coupling rate*)
κR = κS*((c0*ΩR)/ω0)                        (*coupling rate*)

(*SET UP PDEs TO SOLVE*)

Clear[sP, sSt, eqs, vars, bcs];
eqs = {{D[sP[z, t], z] == I κS*(I κR Exp[-γR*t/2]*
Integrate[sP[z, d]*Conjugate[sSt[z, d]]*Exp[γR*d/2], {d, 0, t}])*sSt[z, t]},
{D[sSt[z, t], z] == I κS*Conjugate[I κR Exp[-γR*t/2]*
Integrate[sP[z, d]*Conjugate[sSt[z, d]]*Exp[γR*d/2], {d, 0, t}]]*sP[z, t]}}
vars = {sP[z, t], sSt[z, t]};      (*variables to solve for*)
bcs = {sP[0, t] == Sqrt[0.9], sSt[0, t] == Sqrt[0.1]};   (*boundary conditions at z=0*);
NDSolve[{eqs, bcs}, vars, {t, 0, 0.02}, {z, 0, 0.2}]


EDIT 2: Defining the integration of R as a function outside the NDsolve-block also seems to work.

(*DEFINE ALL THE CONSTANTS*)

Clear[ΩR, αS, κR, κS, γR, sar, c0];
ω0 = 2*Pi*2.99792458*10^8/(1.55*10^(-6));   (*optical frequency in Hz*)
c0 = 2.99792458*10^8;                       (*speed of light in vacuum*)
ΩR = 2*Pi*5.611*10^6;                       (*mechanical frequency in Hz*)
γR = 2*Pi*400;                              (*damping rate/linewidth in Hz*)
κS = -219.666*0.1;                          (*coupling rate*)
κR = κS*((c0*ΩR)/ω0)                        (*coupling rate*)

R[t_, P_, S_] := I κR Exp[-γR*t/2]* Integrate[P*Conjugate[S]*Exp[γR*d/2], {d, 0, t}]

(*SET UP PDEs TO SOLVE*)

Clear[sP, sSt, eqs, vars, bcs];
eqs = {{D[sP[z, t], z] ==
I κS*(R[t, sP[z, t], sSt[z, t]])*sSt[z, t]}, {D[sSt[z, t],
z] == I κS*Conjugate[R[t, sP[z, t], sSt[z, t]]]*sP[z, t]}}
vars = {sP[z, t],sSt[z, t]};(*variables to solve for*)
bcs = {sP[0, t] == Sqrt[0.9], sSt[0, t] == Sqrt[0.1]};(*boundary conditions for input fields at z=0*);
sol = NDSolve[{eqs, bcs}, vars, {t, 0, 0.003}, {z, 0, 0.1}]

(* GENERATE PLOTS *)

DensityPlot[sP[z, t]^2 /. sol, {t, 0, 0.003}, {z, 0, 0.1},
ColorFunction -> "Rainbow",
FrameLabel -> {"Time (s)", "Propagation distance (m)"},
PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 80]

DensityPlot[sSt[z, t]^2 /. sol, {t, 0, 0.003}, {z, 0, 0.1},
ColorFunction -> "Rainbow",
FrameLabel -> {"Time (s)", "Propagation distance (m)"},
PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 80]

DensityPlot[
Abs[R[t, sP[z, t], sSt[z, t]]] /. sol, {t, 0, 0.003}, {z, 0, 0.1},
ColorFunction -> "Rainbow",
FrameLabel -> {"Time (s)", "Propagation distance (m)"},
PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 40]

• Welcome to Mma SE Make the most of the site and take the tour now. Help us to help you, write an excellent question. Edit if improvable, show due diligence, give brief context, include minimal working examples of code and data in formatted form. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. – rhermans May 16 '17 at 14:26
• Here its considered helpful and polite show you own efforts and share your data and code attempts in a well formatted form, so we can quickly see the problem you are facing. Please help us to help you and edit your question to add the code of your own attempts so we can reproduce the error you mentioned. – rhermans May 16 '17 at 14:27
• Your solution is a clever approach but may not be accurate everywhere in {z, t}. This is because, as NDSolve calls the function to compute R, it passes values of var that have been updated for some t and values that have not been updated for others. (Which is which depends on the details of the internal workings of NDSiolve.) I compare our two answers and found differences of as large as about 10%. – bbgodfrey May 17 '17 at 14:06

Although NDSolve cannot solve this system of equations directly, it can produce the desired answer with some help. Here, we take advantage of the small size of the coupling to solve the system iteratively. First, initialize the problem.

ω0 = 2*Pi*2.99792458*10^8/(1.55*10^(-6));   (*optical frequency in Hz*)
c0 = 2.99792458*10^8;                   (*speed of light in vacuum*)
ΩR = 2*Pi*5.611*10^6;     (*mechanical frequency in Hz*)
γR = 2*Pi*400;                   (*damping rate/linewidth in Hz*)
κS = -219.666*0.1;               (*coupling rate*)
κR = κS*((c0*ΩR)/ω0);     (*coupling rate*)

eqsz = {{D[sP[z, t], z] == 0}, {D[sSt[z, t], z] == 0}};
vars = {sP[z, t], sSt[z, t]};
bcs = {sP[0, t] == Sqrt[0.9], sSt[0, t] == Sqrt[0.1]};
sz = vars /. Flatten@NDSolve[{eqsz, bcs}, vars, {t, 0, 0.02}, {z, 0, 0.2}];


Then solve the differential equations alternately for r[z, t] and for {sP[z, t], sSt[z, t]}. Six iterations is more than adequate for good convergence.

Do[
eqst = D[r[z, t], t] + r[z, t] γR/2 == I κR First[sz] Conjugate[Last[sz]];
sr = r[z, t] /. Flatten@NDSolve[{eqst, r[z, 0] == 0}, r[z, t], {t, 0, 0.02}, {z, 0, 0.2}];

eqsz = {{D[sP[z, t], z] == I κS*sr*sSt[z, t]},
{D[sSt[z, t], z] == I κS*Conjugate[sr]*sP[z, t]}};
bcs = {sP[0, t] == Sqrt[0.9], sSt[0, t] == Sqrt[0.1]};
sz = vars /. Flatten@NDSolve[{eqsz, bcs}, vars, {t, 0, 0.02}, {z, 0, 0.2}];,
{i, 6}]

Plot3D[Im@sr, {t, 0, 0.02}, {z, 0, 0.2}, PlotLabel -> R,
AxesLabel -> {t, z}, PlotRange -> All]
Plot3D[sz, {t, 0, 0.02}, {z, 0, 0.2}, PlotLabel -> "sP, sSt",
AxesLabel -> {t, z}]


The computation takes a few seconds.

• Thank you very much! This works very well! I also tried another idea: I've simply defined the integration of R as a function outside the equations used for NDSolve (see attached code in my initial post). – gersch07 May 17 '17 at 11:04