How to solve this trigonometric complex ODE system?

The system of nonlinear ODE is $$\mathrm{i}\,s(\ddot p-\frac{1}{2}\sin{2p}\;\dot q^2)=\mathrm{i}\,c\sin p\;\dot q-a\sin p+b\cos p\cos q\,,\\ \mathrm{i}\,s(\sin^2p\;\ddot q+\sin{2p}\;\dot p\dot q)=-\mathrm{i}\,c\sin p\;\dot p-b\sin{p}\sin{q},$$ wherein $t\in[0,T],T>0$. Also, we have $a,b>0,c=1/2$. The positive $s$ must be small enough since the problem mathematically cares about the $s\rightarrow0$ limit. One can regard it as some strange motion on a sphere and $p,q$ as the two spherical coordinate angles. (So it would be very interesting to visualize any solution on the sphere.) Note that they can in general be complex (hence I said strange) because of the presence of $\mathrm{i}$, although set to be real at the boundaries.

One can scale the equations to have $a^2+b^2=1$. Then the interested domain of parameters is $c=1/2$, $s\ll1=(a^2+b^2)^{-1/2}\ll T$. For instance, let's take $T=6\sim20,s=0.03\sim0.08$ for a mild understanding of $\ll$ and the small $s$ limit. One can tune, of course.

The solution is to some extent expected to exhibit three regions in $t$. Near the left boundary $t-0 \sim O(s)$, rapidly varying from $(p,q)_{t=0}$ to some unknown $(\bar p_1,\bar q_1)$. Near the right boundary $T-t \sim O(s)$, rapidly varying from some unknown $(\bar p_2,\bar q_2)$ to $(p,q)_{t=T}$. And intermediately, slowly varying from $(\bar p_1,\bar q_2)$ to $(\bar p_2,\bar q_2)$. It's possible that $(\bar p_1,\bar q_1)=(\bar p_2,\bar q_2)$ in some cases. Usually, there are more than one valid solutions.

For instance, we hope to solve the following two representative cases $$a=1,b=0,p(0)=p(T)=2\pi/3,q(0)=q(T)=\pi/4,\\ a=1/2,b=\sqrt{3}/2,p(0)=p(T)=3\pi/4,q(0)=\pi/4,q(T)=3\pi/4.$$ I tried NDSolve as follows without any success.

max = 10; s = 0.06; a = 1; b = 0; c = 0.5;
s = NDSolve[{I s (p''[t] - 1/2 Sin[2 p[t]] (q'[t])^2) ==
I c Sin[p[t]] q'[t] - a Sin[p[t]] + b Cos[p[t]] Cos[q[t]],
I s (Sin[p[t]]^2 q''[t] + Sin[2 p[t]] p'[t] q'[t]) == -I c Sin[
p[t]] p'[t] - b Sin[p[t]] Sin[q[t]], p == 2 \[Pi]/3,
p[max] == 2 \[Pi]/3, q == \[Pi]/4, q[max] == \[Pi]/4}, {p,
q}, {t, 0, max}]


A more demanding question

Ultimately, I hope to solve the case with following additional integral terms (setting, e.g., all $K_i=0.1,0.01$) added to the right hand side of the above two equations, respectively. It briefly damps the above motion. Either the case with or without integrals is interested in. If you can solve both, certainly wonderful. But solving only the above pure ODE is good enough for posting an answer. I just feel that it mightn't be good to raise a separate question for the integral case. $$\int_{0}^T{\frac{\mathrm{d}t'}{(t-t')^2}[K_1\cos{p_t} \cos{q_{t}}\sin{p_{t'}} \cos{q_{t'}}+K_2\cos{p_t} \sin{q_{t}}\sin{p_{t'}} \sin{q_{t'}}-K_3\sin{p_t} \cos{p_{t'}}]}\\ \int_{0}^T{\frac{\mathrm{d}t'}{(t-t')^2}[-K_1\sin{p_t} \sin{q_{t}}\sin{p_{t'}} \cos{q_{t'}}+K_2\sin{p_t} \cos{q_{t}}\sin{p_{t'}} \sin{q_{t'}}]}$$

• Maybe you can add some background information so we can analyse your problem easier. – xzczd Dec 10 '17 at 5:38
• @xzczd Please see my updates. BTW, I tried the finite difference method learnt from you. (Thanks again!) It probably works to some extent. But the initial change at any boundary is very very fast and abrupt. Not sure good or not. Perhaps a nonuniform grid (dense around boundary) is better? It is certainly desirable to have other solutions or the FDM solution from much more experienced practitioners like you. In any case, this is a very solid and interesting physical problem. I think it might be good to raise here. – xiaohuamao Dec 10 '17 at 6:50
• There is a dismatch with the initial conditions. Regard the real part of the two equations and ComplexExpand it to get -b Cos[p[t]] Cos[q[t]] + a Sin[p[t]]==0and b Sin[p[t]] Sin[q[t]]==0. Solve this for the intialt =0 p[t]->2 Pi/3 and q[t]->Pi/4to get {{a -> 0, b -> 0}} – Akku14 Dec 10 '17 at 10:19
• @Akku14 No, there's no mismatch. Note that derivatives are not necessarily real at the boundaries. – xiaohuamao Dec 10 '17 at 10:46
• Translating it into Maple and making the substituion $p=u+iv, q = \xi +i \eta$ and putting $,\max =m,m = 0.5$ instead of $\max=10$ (The $\max$ name is a protected name in Maple.), I succeed to solve the problem numerically. See dropbox.com/s/terythorjg7s567/NODES.pdf?dl=0 for details. – user64494 Dec 10 '17 at 18:14

Again, let me give a solution based on finite difference method (FDM). I'll use pdetoae for the generation of difference equation.

Probably you've already noticed that, finding proper initial guess for this problem is no longer as easy as in e.g. your previous problem. We need to guess it in a cleverer way. Inspired by the observations of user64494 and bbgodfrey i.e. the system can be solved more easily when max isn't that large, I found that using solutions at smaller max as initial guess seems to work wonders.

Clear@max
{s = 1/100, c = 1/2, points = 900, difforder = 4};
initialguess = ConstantArray[1/2, points];
{a, b, {{p0, pmax}, {q0, qmax}}} = {1, 0, {{2 Pi/3, 2 Pi/3}, {Pi/4, Pi/4}}};
{eq, bc} = {{I s (p''[t] - 1/2 Sin[2 p[t]] (q'[t])^2) ==
I c Sin[p[t]] q'[t] - a Sin[p[t]] + b Cos[p[t]] Cos[q[t]],
I s (Sin[p[t]]^2 q''[t] + Sin[2 p[t]] p'[t] q'[t]) == -I c Sin[p[t]] p'[t] -
b Sin[p[t]] Sin[q[t]]}, {p == p0, p[max] == pmax, q == q0, q[max] == qmax}};

grid = Array[# &, points, {0, max}];
ptoafunc = pdetoae[{p, q}[t], grid, difforder];
delete = #[[2 ;; -2]] &;
ae = delete /@ ptoafunc[eq];
solandmaxlst = NestList[Function[guessandmax, Block[{guess, max, guessfunc, sollst},
{guess, max} = guessandmax;
guessfunc = ListInterpolation[#, {0, max}] & /@ guess;
sollst =
ArrayReshape[
FindRoot[{ae, bc},
Flatten[#, 1] &@
Table[{{p, q}[[var]][index], guessfunc[[var]][index]}, {var, 2}, {index,
grid}], MaxIterations -> 1000][[All, -1]], {2, points}];
{sollst, max + 1}]], {{initialguess, initialguess}, 1}, 10]; // AbsoluteTiming
(* {21.842398, Null} *)

funclst = Table[
ListInterpolation[#, {0, lst[] - 1}] & /@ lst[], {lst, solandmaxlst // Rest}];


We can observe how the solution evolutes by:

plotlst = Table[
Plot[{Re@#[t], Im@#[t]}, {t, 0, #["Domain"][[1, -1]]}, PlotRange -> All] & /@ pq, {pq,
funclst}];
ListAnimate@plotlst The final solution is

plotlst[[-1]] // GraphicsColumn Further check by adjusting points etc. shows the solution seems to be reliable.

Finally the following is the solution for {a, b, {{p0, pmax}, {q0, qmax}}} = {1/2, Sqrt/2, {{3 Pi/4, 3 Pi/4}, {Pi/4, 3 Pi/4}}}: • Thanks for your answer! I clarified the parameters a bit. Perhaps it helps. – xiaohuamao Dec 12 '17 at 9:37
• Learned again the importance of good initial guess. BTW, I previously mentioned trying to use the Logistic function to generate a nonuniform grid (sparse in middle) for the FDM on the pure ODE. But it seems that even if a tiny nonuniformity slows down the solution a lot, which doesn't happen in your previous Gauss rule grid. Do you know any reason or have any comments. Like this one. Logistic[x_, s_] = 1/(1 + Exp[-(x - 1)/(10 Sqrt[s])]); grid = Rescale[Table[Logistic[x, s], {x, 0, 2, 2/points}], {Logistic[0, s], Logistic[2, s]}, domain]; – xiaohuamao Dec 12 '17 at 10:24
• +1. It's me. The question arises: why isn't it implemented in Mathematica yet? – user64494 Dec 12 '17 at 10:48
• @xiaohuamiao I'm not sure if it's the whole reason, but it's related to the issue mentioned in this post. Specifically speaking, the exact form of your grid is complicated, so the deduced difference equation becomes complicated and hard for Mathematica to handle. If you modify the code to SetAttributes[#, NHoldAll] & /@ {p, q}; ……FindRoot[N@{ae, bc}…… it'll be much faster. – xzczd Dec 12 '17 at 12:39
• @user64494 I sincerely believe it's because of some deep reasons that are beyond the understanding of us average users, rather than something like they think FDM is… not cool. – xzczd Dec 12 '17 at 12:46

It is convenient and also provides additional insight to rescale t by s, for which the ODE system becomes,

ss = NDSolveValue[{(p''[t] - 1/2 Sin[2 p[t]] (q'[t])^2) == c Sin[p[t]] q'[t] -
a s Sin[p[t]]/I,
(Sin[p[t]]^2 q''[t] + Sin[2 p[t]] p'[t] q'[t]) == -c Sin[p[t]] p'[t],
p == 2 π/3, p[max/s] == 2 π/3, q == π/4, q[max/s] == π/4}, {p, q}, {t, 0, max/s}]
// Flatten;


Thus, s now can be absorbed into a in the first ODE, if desired, and does not appear at all in the second ODE. Note, though, that the range of integration now is {0, max/s}, so s is not completely eliminated. (b has been dropped here, because it is set to zero in the question. However, it can be reintroduced easily, if desired.) It also is important to observe that the coefficient of q''[t] is proportional to Sin[p[t]]^2. Thus, the ODE system becomes singular whenever p becomes equal to π anywhere in the domain of integration.

User64494 used Maple to obtain a solution for

max = 1/2; s = 1/100; a = 1; b = 0; c = 1/2;


It can be obtained with Mathematica by

ss = NDSolveValue[{(p''[t] - 1/2 Sin[2 p[t]] (q'[t])^2) == c Sin[p[t]] q'[t] -
a s Sin[p[t]]/I,
(Sin[p[t]]^2 q''[t] + Sin[2 p[t]] p'[t] q'[t]) == -c Sin[p[t]] p'[t],
p == 2 π/3, p[max/s] == 2 π/3, q == π/4, q[max/s] == π/4}, {p, q},
{t, 0, max/s}, Method -> {"Shooting", "ImplicitSolver" -> {"Newton",
"StepControl" -> "LineSearch"}, "StartingInitialConditions" ->
p' == 1/100 - 20/100 I, q' == 1/100}}] // Flatten;
Plot[Evaluate@ReIm@Through[ss[t]], {t, 0, max/s}] Here, the real and imaginary parts of p are blue and orange, and of q are green and red. (See the answer to 124457 for a discussion of "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}, without which NDSolve crashes.) With modest effort max can be increased to 55/100 by means of "StartingInitialConditions" -> {p' == 4/100 - 29/100 I, q' == 1/100 - 8/100 I}. Visibly, the ODE system almost is singular, and I have been unable to obtain a solution for max == 56/100. There may, of course, be solutions for stable islands of larger max, but I have not tried to find them. In that regard, it is easy to obtain solutions for max == 40/100 and max == 42/100 but not for max == 41/100, again because the ODE system is singular there.

An alternative approach to finding a solution for the specific parameters in the question is to begin with small a and gradually increase it. Unfortunately, I was unable to increase a beyond 4 10^-2 (and max == 10), for which with "StartingInitialConditions" -> {p' == 0.00026069739639169307 - 0.1478258792803156 I, q' == -0.1646363401819224 + 0.01211739772818876 I}, the solution is Larger values of a result in

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

but increasing WorkingPrecision` does not help and is painfully slow besides. In closing, I have been able to find solutions for many parameters but not for those in the question.

• Thank you a lot for your nice and detailed answer. The only crucial problem is max less than 1 is way smaller than that needed to show the expected features. And I actually already rescaled the equations to make $a^2+b^2=1$, so the 'increasing $a$' approach is somewhat not very good. But I learned and enjoyed your answer. – xiaohuamao Dec 12 '17 at 5:52
• I clarified the parameters a bit. Perhaps it helps. – xiaohuamao Dec 12 '17 at 9:38