help with complex second order differential equation and NDSolve

I have a complicated 2nd ODE with boundary values. It is actually a vector so there are three ODEs and they are coupled. I would like to solve for sx[x], sy[x], sz[x] with the boundaries sx[0] = 1 and sy[0] = 0, and sz[0] = 0. Also alsosx[L] = sy[L] = sz[L] = 0 where L should be infinity but I'm not sure if that's allowed. So I just take it to be a large number. Note that the s[x]'s exist in the denominator of one of the terms and are non-linear.

I have tried this with NDSolve as you can see below. For some values of the parameters I can achieve results in reasonable amount of time. But if I make, for instance, Bz too large (like 10), mathematica slows down considerably and I may get errors.

I would like to know if I can hope to use NDSolve to find solutions for a broad range of input parameters. And if not, I suspect that a finite different method must be used. I have seen a canned finite difference method for ODEs but it was not general enough to account for the complex form and nonlinearity that I have here. If someone could help me get started on putting this into a finite difference framework, or pointing me to good resources, that would also be appreciated.

B = {0, 0, Bz};

params = {Ex -> 1., m -> 0.5, d -> 4.,
Bz -> 1., g -> 1., tau -> 1., be -> 0.1, BL -> 0.01,
bN -> -1, L -> 20., s0 -> 1.};

sol = Flatten[
NDSolve[{d {sx''[x], sy''[x], sz''[x]} + m Ex {sx'[x],
sy'[x], sz'[x]} + g  Cross[
B + bN  #.(B + be #)/((B + be #).(B + be #) +
BL^2) (B + be #), #] - #/tau == {0., 0., 0.},
sx[L] == 0., sy[L] == 0., sz[L] == 0., sx[0.] == s0,
sy[0.] == 0., sz[0.] == 0.} /. params, {sx[x], sy[x],
sz[x]}, {x, 0, 20.}] ] &[{sx[x], sy[x], sz[x]}];

Sx = sx[x] /. sol;
Sy = sy[x] /. sol;
Sz = sz[x] /. sol;

• Cross[B + bN #.(B + be #)/((B + be #).(B + be #) + BL^2) (B + be #), #] &[{sx[x], sy[x], sz[x]}] has a LeafCount of 3843. After Simplify it drops to 191, suggesting that the ODEs should be simplified before being inserted into NDSolve. Also, for the parameters given sz[x] is zero. Deleting it from the set of ODEs could improve speed somewhat. – bbgodfrey Jan 6 '17 at 0:46
• The result of these changes is to reduce runtime from 146 sec to 18 sec. – bbgodfrey Jan 6 '17 at 1:12
• sz[x] will not be zero in general. But thanks for the LeafCount info - was not aware of that function. – BeauGeste Jan 6 '17 at 1:30
• @bbgodfrey When I run this I am getting an error in apparently in FindRoot that it is maxed out at 100 iterations. I am not using FindRoot but I presume that NDSolve is. How can I increase the iterations in this case? – BeauGeste Jan 6 '17 at 1:35
• According to my limited experience, finite difference method (FDM) is more likely to fail when dealing with boundary value problem (BVP) of nonlinear ODE. (It needs a good enough initial guess in the whole domain! ) If your problem can't be handled by NDSolve, I'm afraid FDM won't help, either. (FDM will always give you an answer of course, but it's probably wrong. ) – xzczd Jan 6 '17 at 3:02

Numerical Solution

The following actions greatly improve speed and accuracy for the particular ODE system in the question. First, the function involving

Cross[B + bN  #.(B + be #)/((B + be #).(B + be #) + BL^2) (B + be #), #] &
[{sx[x], sy[x], sz[x]}]


should be performed outside NDSolve and Simplify applied to the result to reduce running time by an order of magnitude or so. Moreover, with B perpendicular to Ex, as is the case here, the component of s parallel to B, here sz, is zero and can be omitted from the problem. Next, use the "Shooting" Method with increased "MaxIterations". Finally, as Bz is increased, L can be decreased. With these changes, the problem with Bz == 10 becomes

sol = Flatten@NDSolve[{d {sx''[x], sy''[x]} + m Ex {sx'[x], sy'[x]} +
g  {-((Bz sy[x] (BL^2 + Bz^2 + be (be + bN) sx[x]^2 +
be (be + bN) sy[x]^2))/(BL^2 + Bz^2 + be^2 sx[x]^2 + be^2 sy[x]^2)),
(Bz sx[x] (BL^2 + Bz^2 + be (be + bN) sx[x]^2 +
be (be + bN) sy[x]^2))/(BL^2 + Bz^2 + be^2 sx[x]^2 + be^2 sy[x]^2)} -
{sx[x], sy[x]}/tau == {0, 0},
sx[L] == 0, sy[L] == 0, sx[0] == s0, sy[0] == 0} /. params,
{sx[x], sy[x]}, {x, 0, 10},
Method -> {"Shooting", "StartingInitialConditions" -> {sx[0] == s0 /. params,
sy[0] == 0, sx'[0] == -1, sy'[0] == 1}, "MaxIterations" -> 500}];


and produces a result in seconds.

Plot[{sx[x], sy[x]} /. sol, {x, 0, 10}, PlotRange -> All]


Approximate Symbolic Solution

It is worth noting that the ODE system can be solved approximately with DSolve for large Bz or for large x.

asol = Flatten@DSolve[{d {sx''[x], sy''[x]} + m Ex {sx'[x], sy'[x]} +
g {-Bz sy[x] , Bz sx[x]} - {sx[x], sy[x]}/tau == {0, 0}}, {sx[x], sy[x]}, x] // Simplify;


Although the result is far too lengthy to display here, its exponential factors can be isolated without difficulty.

aexp = Simplify[Union@Cases[asol // Expand, Exp[_], Infinity]]
(* {E^(-((Ex m x + (Sqrt[4 d + 4 I Bz d g tau + Ex^2 m^2 tau] x)/Sqrt[tau])/(2 d))),
E^(-((Ex m x + (Sqrt[4 d - 4 I Bz d g tau + Ex^2 m^2 tau] x)/Sqrt[tau])/(2 d))),
E^(-((Ex m x - (Sqrt[4 d + 4 I Bz d g tau + Ex^2 m^2 tau] x)/Sqrt[tau])/(2 d))),
E^(-((Ex m x - (Sqrt[4 d - 4 I Bz d g tau + Ex^2 m^2 tau] x)/Sqrt[tau])/(2 d)))} *)


Only the exponentially decaying terms satisfy the boundary conditions at large L. The growing solutions can be eliminated as follows.

asolsx = Collect[((sx[x] /. asol) // Expand) /.
Exp[z_] :> Exp[Simplify[z]], {Exp[_], C[_]}, Simplify];
asolsy = Collect[((sy[x] /. asol) // Expand) /.
Exp[z_] :> Exp[Simplify[z]], {Exp[_], C[_]}, Simplify];
c24 = Collect[Flatten@Solve[{Coefficient[asolsx, aexp[[3]]] == 0,
Coefficient[asolsx, aexp[[4]]] == 0}, {C[2], C[4]}], C[_], Simplify];
asolsx1 = Collect[asolsx /. c24, {Exp[_], C[_]}, Simplify]
asolsy1 = Collect[asolsy /. c24, {Exp[_], C[_]}, Simplify]
(* aexp[[1]] (C[1]/2 - 1/2 I C[3]) + aexp[[2]] (C[1]/2 + 1/2 I C[3]) *)
(* aexp[[2]] (-(1/2) I C[1] + C[3]/2) + aexp[[1]] (1/2 I C[1] + C[3]/2) *)


where aexp[[1]] and aexp[[2]] have been written in place of the corresponding exponentials solely for notational simplicity. The remaining constants are evaluated from

{asolsx1 == 1, asolsy1 == 0} /. x -> 0 /. Equal -> Rule
(* {C[1] -> 1, C[3] -> 0} *)


and

Plot[{asolsx1, asolsy1} /. % /. params, {x, 0, 10}, PlotRange -> All]


gives a plot indistinguishable to the eye from the plot above. Thus, we have a good symbolic solution for BL^2 + Bz^2 large compared be^2 and Abs[be (be + bN)].

When this is not the case, these solutions instead can be used to obtain superior boundary conditions at large L, allowing L to be reduced in magnitude. Doing so is useful, because numerical solutions with NDSolve inevitably pick up some of the exponentially growing solutions, causing loss of precision in fitting the boundary conditions at L. To generate these boundary conditions, first construct

Simplify[asolsx1 + I asolsy1]
(* aexp[[2]] (C[1] + I C[3]) *)


and take its logarithmic derivative.

bc = D[Log[sx[x] + I sy[x]], x] == D[Log[%], x]
(* (Derivative[1][sx][x] + I Derivative[1][sy][x])/(sx[x] + I sy[x]) ==
-((Ex m + Sqrt[4 d - 4 I Bz d g tau + Ex^2 m^2 tau]/Sqrt[tau])/(2 d)) *)


The Real and Imaginary parts of bc are the desired boundary conditions. For instance, with Bz == 10, as above,

sx[x] + I sy[x] == ComplexExpand[(Derivative[1][sx][x] + I Derivative[1][sy][x])/
N[Last@bc /. params]]
(* sx[x] + I sy[x] == -0.674697 Derivative[1][sx][x] + I (-0.514645 Derivative[1][sx][x] -
0.674697 Derivative[1][sy][x]) + 0.514645 Derivative[1][sy][x] *)
{bc1, bc2} = Thread[{sx[x], sy[x]} == Refine[ReIm[%[[2]]],
Derivative[1][sx][x] | Derivative[1][sy][x]) ∈ Reals]] /. x -> L
(* {sx[L] == -0.674697 Derivative[1][sx][L] + 0.514645 Derivative[1][sy][L],
sy[L] == -0.514645 Derivative[1][sx][L] - 0.674697 Derivative[1][sy][L]} *)