1
$\begingroup$

I have four coupled first-order non-linear differential equations, denoted as: A1[x], A2[x], A3[x], A4[x] which are all functions of x. I have the following code which attempts to solve the equations using ParametricSolveND by varying one of the initial conditions of the parameter (namely, A4[0] which I have denoted as the parameter j).

ω1 = 2 π*5*10^9;
ω2 = 2 π*5*10^9;
ω3 = 2 π*3*10^9;
ω4 = ω1 + ω2 - ω3;
Cj = 329*10^-15;
LL = 100*10^-12;
a = 10*10^-6;
I0 = 3.29*10^-6;
CC0 = 39*10^-15;
k1 = (Sqrt[CC0 LL] *(ω1))/(a Sqrt[1 - Cj LL *(ω1)^2]);
k2 = (Sqrt[CC0 LL] *(ω2))/(a Sqrt[1 - Cj LL *(ω2)^2]);
k3 = (Sqrt[CC0 LL] *(ω3))/(a Sqrt[1 - Cj LL *(ω3)^2]);
k4 = (Sqrt[CC0 LL] *(ω4))/(a Sqrt[1 - Cj LL *(ω4)^2]);
Δkl = k1 + k2 - k3 - k4;
κ1 = (a^4*k1*k2*k3*k4*(k3 + k4 - k2))/(8*CC0*I0^2*LL^3*ω1^2);
κ2 = (a^4*k1*k2*k3*k4*(k3 + k4 - k1))/(8*CC0*I0^2*LL^3*ω2^2);
κ3 = (a^4*k1*k2*k3*k4*(k1 + k2 - k4))/(8*CC0*I0^2*LL^3*ω3^2);
κ4 = (a^4*k1*k2*k3*k4*(k1 + k2 - k3))/(8*CC0*I0^2*LL^3*ω4^2);
α11 = (a^4*k1^5)/(16*CC0*I0^2*LL^3*ω1^2);
α12 = (a^4*k1^3*k2^2)/(8*CC0*I0^2*LL^3*ω1^2);
α13 = (a^4*k1^3*k3^2)/(8*CC0*I0^2*LL^3*ω1^2);
α14 = (a^4*k1^3*k4^2)/(8*CC0*I0^2*LL^3*ω1^2);
α21 = (a^4*k2^3*k1^2)/(8*CC0*I0^2*LL^3*ω2^2);
α22 = (a^4*k2^5)/(16*CC0*I0^2*LL^3*ω2^2);
α23 = (a^4*k2^3*k3^2)/(8*CC0*I0^2*LL^3*ω2^2);
α24 = (a^4*k2^3*k4^2)/(8*CC0*I0^2*LL^3*ω2^2);
α31 = (a^4*k3^3*k1^2)/(8*CC0*I0^2*LL^3*ω3^2);
α32 = (a^4*k3^3*k2^2)/(8*CC0*I0^2*LL^3*ω3^2);
α33 = (a^4*k3^5)/(16*CC0*I0^2*LL^3*ω3^2);
α34 = (a^4*k3^3*k4^2)/(8*CC0*I0^2*LL^3*ω3^2);
α41 = (a^4*k4^3*k1^2)/(8*CC0*I0^2*LL^3*ω4^2);
α42 = (a^4*k4^3*k2^2)/(8*CC0*I0^2*LL^3*ω4^2);
α43 = (a^4*k4^3*k3^2)/(8*CC0*I0^2*LL^3*ω4^2);
α44 = (a^4*k4^5)/(16*CC0*I0^2*LL^3*ω4^2) // N;

system = {A1'[x] == I*κ1*Conjugate[A2[x]]*A3[x]*A4[x]*E^(-I*Δkl*x) + I*A1[x]*(α11*Abs[A1[x]]^2 + α12*Abs[A2[x]]^2 + α13*Abs[A3[x]]^2 + α14*Abs[A4[x]]^2), 
          A2'[x] == I*κ2*Conjugate[A1[x]]*A3[x]*A4[x]*E^(-I*Δkl*x) + I*A2[x]*(α21*Abs[A1[x]]^2 + α22*Abs[A2[x]]^2 + α23*Abs[A3[x]]^2 + α24*Abs[A4[x]]^2), 
          A3'[x] == I*κ3*A1[x]*A2[x]*Conjugate[A4[x]]*E^(I*Δkl*x) + I*A3[x]*(α31*Abs[A1[x]]^2 + α32*Abs[A2[x]]^2 + α33*Abs[A3[x]]^2 + α34*Abs[A4[x]]^2), 
          A4'[x] == I*κ4*A1[x]*A2[x]*Conjugate[A3[x]]*E^(I*Δkl*x) + I*A4[x]*(α41*Abs[A1[x]]^2 + α42*Abs[A2[x]]^2 + α43*Abs[A3[x]]^2 + α44*Abs[A4[x]]^2), 
          A1[0] == (I0*25)/ω1, A2[0] == (I0*25)/ω2, A3[0] == 0, A4[0] == j};

DEsols = ParametricNDSolve[system, {A1[x], A2[x], A3[x], A4[x]}, {x, 0, 2000}, {j}]
Plot[Evaluate@Table[Abs[(A4[j][x]) /. DEsols]^2, {j, 0, 10}], {x, 0, 2000}, PlotStyle -> {Orange}, PlotLegends -> {"A4"}, PlotRange -> All, AxesOrigin -> {0, 0}]

However, it is not plotting and I'm not sure what I've done wrong. Furthermore, I intend to plot A4[x] as a function of A4[0] for a fixed x (x=2000). How should I go about fixing this? Thank you.

$\endgroup$
1
$\begingroup$

After modifying ParametricNDSolve to ParametricNDSolveValue (without argument brackets ...[x])

DEsols = ParametricNDSolveValue[system, {A1 , A2 , A3 , A4 }, {x, 0, 2000}, {j}]

you can access A4

Plot[sols[0][[4]][x], {x, 0, 2000}]

Unfortunately your system seems to be solvable only for j==0

| improve this answer | |
$\endgroup$
1
$\begingroup$

We can use NDSolve[] to solve this problem as follows:

system[j_] := {A1'[x] == 
    I*\[Kappa]1*Conjugate[A2[x]]*A3[x]*A4[x]*
      E^(-I*\[CapitalDelta]kl*x) + 
     I*A1[x]*(\[Alpha]11*Abs[A1[x]]^2 + \[Alpha]12*
         Abs[A2[x]]^2 + \[Alpha]13*Abs[A3[x]]^2 + \[Alpha]14*
         Abs[A4[x]]^2), 
   A2'[x] == 
    I*\[Kappa]2*Conjugate[A1[x]]*A3[x]*A4[x]*
      E^(-I*\[CapitalDelta]kl*x) + 
     I*A2[x]*(\[Alpha]21*Abs[A1[x]]^2 + \[Alpha]22*
         Abs[A2[x]]^2 + \[Alpha]23*Abs[A3[x]]^2 + \[Alpha]24*
         Abs[A4[x]]^2), 
   A3'[x] == 
    I*\[Kappa]3*A1[x]*A2[x]*Conjugate[A4[x]]*
      E^(I*\[CapitalDelta]kl*x) + 
     I*A3[x]*(\[Alpha]31*Abs[A1[x]]^2 + \[Alpha]32*
         Abs[A2[x]]^2 + \[Alpha]33*Abs[A3[x]]^2 + \[Alpha]34*
         Abs[A4[x]]^2), 
   A4'[x] == 
    I*\[Kappa]4*A1[x]*A2[x]*Conjugate[A3[x]]*
      E^(I*\[CapitalDelta]kl*x) + 
     I*A4[x]*(\[Alpha]41*Abs[A1[x]]^2 + \[Alpha]42*
         Abs[A2[x]]^2 + \[Alpha]43*Abs[A3[x]]^2 + \[Alpha]44*
         Abs[A4[x]]^2), A1[0] == (I0*25)/\[Omega]1, 
   A2[0] == (I0*25)/\[Omega]2, A3[0] == 0, A4[0] == j}; 
sol[j_] := NDSolve[system[j], {A1, A2, A3, A4}, {x, 0, 2000}]

Now we can plot solution for small range of j

 {LogLogPlot[
  Evaluate@Table[
    Abs[(A1[x]) /. sol[j]]^2, {j, .5 10^-10, 
     2 10^-10, .5 10^-10}], {x, 10^-5, 2000}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "A1"], 
 LogLogPlot[
  Evaluate@Table[
    Abs[(A2[x]) /. sol[j]]^2, {j, .5 10^-10, 
     2 10^-10, .5 10^-10}], {x, 10^-5, 2000}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "A2"], 
 LogLogPlot[
  Evaluate@Table[
    Abs[(A3[x]) /. sol[j]]^2, {j, .5 10^-10, 
     2 10^-10, .5 10^-10}], {x, 10^-5, 2000}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "A3"], 
 LogLogPlot[
  Evaluate@Table[
    Abs[(A4[x]) /. sol[j]]^2, {j, .5 10^-10, 
     2 10^-10, .5 10^-10}], {x, 10^-5, 2000}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "A4"]}

Figure 1

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for the response. However when I try to plot A4[2000] as a function of j, a bunch of warnings and errors showed up. In particular, I did LogLinearPlot[Abs[(A4[2000]) /. sol[j]]^2, {j, 0, 50}, PlotRange -> All] and it refuses to plot. $\endgroup$ – kowalski Jul 9 at 16:35
  • $\begingroup$ @kowalski Just check that my plots made for {j, .5 10^-10, 2 10^-10, .5 10^-10}. And for this small parameter we have A4 of order 10^29. With increasing j over 3*10^-6 there is no stable solution. $\endgroup$ – Alex Trounev Jul 9 at 16:59
  • $\begingroup$ I can see from the fourth panel A4 that it is indeed in the order of 10^29. However, if I specify my methods, namely sol[j_] := NDSolve[system[j], {A1, A2, A3, A4}, {x, 0, 2000}, Method -> "Extrapolation", StartingStepSize -> 0.01] over j up to 3*10^-6 I do see a plot but I am unsure if it's the right plot since there are many warnings $\endgroup$ – kowalski Jul 9 at 17:22
  • $\begingroup$ @kowalski What actually do you try to describe by this system? $\endgroup$ – Alex Trounev Jul 9 at 18:18
  • $\begingroup$ A1, A2, A3, A4 are amplitudes of a wave in an optic fiber amplifier. A1 and A2 represents the pump amplitudes, A3 represents the idler amplitude and A4 represents the signal amplitude. I am trying to plot gain which is Log10[Abs[A4[x]/A4[0]]^2] as a function of A4[0] which is j to demonstrate saturation of gain as a function of input signal power (Abs[A4[0]]^2) $\endgroup$ – kowalski Jul 9 at 18:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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