0
$\begingroup$

I faced with unusual (at least for me) problems.

I have a three nonlinear equations

$ \dot x = f(x,y,z)\quad \dot y = g(x,y,z)\quad \dot z = h(x,y,z) $

and the following initial conditions $x(0)=1$, $y(0)=0$, $z(4)=1$.

I try to solve it in the Mathematica

L = {D[Ap[z], z] == -gp/2 (Abs[Af[z]]^2 + Abs[Ab[z]]^2) Ap[z] + 
I gpp (Abs[Ap[z]]^2 + 2 Abs[Af[z]]^2 + 2 Abs[Ab[z]]^2) Ap[z] - 
alp/2 Ap[z], 
D[Af[z], z] == -gs/2 Abs[Ap[z]]^2 Af[z] + 
I gss (2 Abs[Ap[z]]^2 + Abs[Af[z]]^2 + 2 Abs[Ab[z]]^2) Af[z] + 
I kk Ab[z] + I db Af[z] - a1s/2 Af[z], -D[Ab[z], z] == 
gs/2 Abs[Ap[z]]^2 Ab[z] + 
I gss (2 Abs[Ap[z]]^2 + 2 Abs[Af[z]]^2 + Abs[Ab[z]]^2) Ab[z] + 
I kk Af[z] + I db Ab[z] - a1s/2 Ab[z]}

vp = 1.; vs = 2; gp = 3; gs = 4; kk = 5; a1s = 0.5; gpp = 0.1; gss = 0.23; alp = 0.79; db = 0.07;

s = NDSolve[Join[L, {Ap[0] == 1, Af[0] == 0, Ab[4] == 1}], {Ap, Af, Ab}, {z, 0, 4}]

But the solution $s$ is not satisfied the condition $Ab[4] =1$. In the Mathematica solution $Ab[4]=-7.27787 + 1.55824 I$ instead of $1$.

How to solve it correctly?

$\endgroup$
2
  • $\begingroup$ Is the solution defined in the complex plane Z or are only complex functions Ap, Af, Ab of real variable z? $\endgroup$ Oct 7, 2018 at 5:38
  • $\begingroup$ $z$ is a real. Functions $Ap, Af, Ab$ are complex. $\endgroup$
    – Peter
    Oct 7, 2018 at 5:47

1 Answer 1

2
$\begingroup$

First, we reduce the system of equations to a real form. Secondly, we will define the initial data using a parametric solution.

L = {-D[Ap[z], z] - gp/2 (Abs[Af[z]]^2 + Abs[Ab[z]]^2) Ap[z] +
     I gpp (Abs[Ap[z]]^2 + 2 Abs[Af[z]]^2 + 2 Abs[Ab[z]]^2) Ap[z] - 
    alp/2 Ap[z], -D[Af[z], z] - gs/2 Abs[Ap[z]]^2 Af[z] + 
    I gss (2 Abs[Ap[z]]^2 + Abs[Af[z]]^2 + 2 Abs[Ab[z]]^2) Af[z] + 
    I kk Ab[z] + I db Af[z] - a1s/2 Af[z], 
   D[Ab[z], z] + gs/2 Abs[Ap[z]]^2 Ab[z] + 
    I gss (2 Abs[Ap[z]]^2 + 2 Abs[Af[z]]^2 + Abs[Ab[z]]^2) Ab[z] + 
    I kk Af[z] + I db Ab[z] - a1s/2 Ab[z]};

Ap'[z] = Ap1'[z] + I*Ap2'[z];

Af'[z] = Af1'[z] + I*Af2'[z];
Ab'[z] = Ab1'[z] + I*Ab2'[z];
 Ap[z] = Ap1[z] + I*Ap2[z];
   Af[z] = Af1[z] + I*Af2[z];
Ab[z] = Ab1[z] + I*Ab2[z];

L // FullSimplify;
eq = ComplexExpand[L];
 eqs = 
  Table[{ComplexExpand[Re[eq[[i]]]] == 0, 
    ComplexExpand[Im[eq[[i]]]] == 0}, {i, 1, 3}];

 vp = 1.; vs = 2; gp = 3; gs = 4; kk = 5; a1s = 0.5; gpp = \
0.1; gss = 0.23; alp = 0.79; db = 0.07;

AB = 
  ParametricNDSolveValue[
   Join[eqs, {Ap1[0] == 1, Af1[0] == 0, Ab1[0] == p, Ap2[0] == 0, 
     Af2[0] == 0, Ab2[0] ==q}], Ab1, {z, 0, 4}, {p,q}];
AC = ParametricNDSolveValue[Join[eqs, {Ap1[0] == 1, Af1[0] == 0, Ab1[0] == 
p, Ap2[0] == 0, Af2[0] == 0, Ab2[0] == q}], Ab2, {z, 0, 4}, {p, q}];
   FindRoot[{AB[p, q][4] == 1, AC[p, q][4] == 0}, {p, -0.5,
0}, {q, -.5, 0}, Method -> "Secant"]

Out[]= {p -> 4.3433*10^-8, q -> 5.71603*10^-9}

 {Plot[Evaluate[{AB[4.3433002424791235`*^-8, 5.716028592104211`*^-9][
 z], AC[4.3433002424791235`*^-8, 5.716028592104211`*^-9][z]}], {z,
0, 4}, PlotRange -> All], 
 Plot[Evaluate[{AB[4.3433002424791235`*^-8, 5.716028592104211`*^-9][
     z], AC[4.3433002424791235`*^-8, 5.716028592104211`*^-9][z]}], {z,
    0, 4}, PlotRange -> Automatic, PlotLegends -> Automatic]}

fig1

Finally, we find the entire solution.

s = NDSolve[
  Join[eqs, {Ap1[0] == 1, Af1[0] == 0, 
    Ab1[0] == 4.3433002424791235`*^-8, Ap2[0] == 0, Af2[0] == 0, 
    Ab2[0] == 5.716028592104211`*^-9}], {Ap1, Ap2, Af1, Af2, Ab1, 
   Ab2}, {z, 0, 4}]
 {Plot[Evaluate[{Ap1[z], Ap2[z], Af1[z], Af2[z], Ab1[z], Ab2[z]} /. 
    s], {z, 0, 4}, PlotRange -> {-.1, 0.3}], 
 Plot[Evaluate[{Ap1[z], Ap2[z], Af1[z], Af2[z], Ab1[z], Ab2[z]} /. 
    s], {z, 0, 4}, PlotLegends -> Automatic, PlotRange -> All]}

fig2

$\endgroup$
7
  • $\begingroup$ Thank you! I understand the idea. But realization of this idea is not right. Because in the function ParametricNDSolveValue should be to parameters $p1$ and $p2$. Ab1[0]=p1 and Ab2[0]=p2. In the functions FindRoot will be two equation $Ab1[4]==1$ and $Ab2[4]==0$. But I faced with the problem that FindRoot does not give me a result and says "FindRoot::nlnum: "The function value {-1.+0.1[4.],0.2[4.]} is not a list of numbers with dimensions {2} at {p1,p2} = {0.1,0.2}."" $\endgroup$
    – Peter
    Oct 7, 2018 at 15:20
  • $\begingroup$ @Peter In your problem, there are no boundary conditions for the imaginary parts of the functions. Think about what the conditions should be, and then I will help solve the problem. $\endgroup$ Oct 7, 2018 at 15:42
  • $\begingroup$ In my problem Ab(4)=1. It means that Ab1(4)=1 and Ab2(4)=0. In your answers Ab2(4) has non zero value. $\endgroup$
    – Peter
    Oct 7, 2018 at 15:52
  • $\begingroup$ Are other boundary conditions correct? $\endgroup$ Oct 7, 2018 at 16:12
  • $\begingroup$ Yes. The boundary conditions are $Ap[0] == 1, \quad Af[0] == 0,\quad Ab[4] == 1$. In the your notations are $Ap1[0] == 1,\quad Af1[0] == 0, \quad Ab1[4] == 1 \quad Ap2[0] == 0, \quad Af2[0] == 0, \quad Ab2[4] == 0$ $\endgroup$
    – Peter
    Oct 7, 2018 at 16:37

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.