3
$\begingroup$

I am trying to solve a system of coupled differential equations

xmax = 15;ϵ = 0.1;
s = NDSolve[{z''[x] == (1 - ϵ/2)*Sin[z[x]] - ϵ/2*Sin[y[x]], 
y''[x] == (1 - ϵ/2)*Sin[y[x]] - ϵ/2*Sin[z[x]], 
y[-xmax] == 0, y[xmax] == 2*π, z[-xmax] == 0,z[xmax] == 0}, {y, z}, {x, -xmax, xmax}, 
MaxSteps -> 10^8, AccuracyGoal -> Automatic, PrecisionGoal -> 50,
WorkingPrecision -> MachinePrecision];
Plot[0.5*(y[x] + z[x]) /. s, {x, -xmax, xmax}, PlotRange -> All] //Print;

If I run this for small value of xmax like xmax = 5, it works just fine and gives the solution as expected. But for higher values of xmax, it runs into tolerance problems and doesn't give the solution as expected. What is the reason for this error?

I see now that for small $\epsilon$, I can uncouple the equation and solve the following equation but it still runs into the same tolerance problem:

 xmax = 5; ϵ = 0;
 s = NDSolve[{y''[x] == y[x] - 0.5*Sin[4*ArcTan[Exp[-x]]],y[-xmax] == 0, y[xmax] == 0}, y, 
{x, -xmax, xmax},MaxSteps -> 10^8, AccuracyGoal -> 8, PrecisionGoal -> 30, 
     WorkingPrecision -> MachinePrecision];
 Plot[y[x] /. s, {x, -xmax, xmax}, PlotRange -> All] // Print;

How can I correct this?

$\endgroup$
2
  • $\begingroup$ Are you actually approximate $\infty$ with xmax? $\endgroup$
    – xzczd
    Jan 26, 2017 at 10:11
  • $\begingroup$ @xzczd Yes, I am approaching infinity with $xmax$. $\endgroup$
    – cleanplay
    Jan 26, 2017 at 13:12

2 Answers 2

3
$\begingroup$

Comment

The solution to the problem BVP system seems to be highly oscillatory. I tried to solve it with Shooting for different combination of StartingInitialConditions. The results are promising but not 100% accurate.

Eq1 = z''[x] == (1 - ϵ/2)*Sin[z[x]] -ϵ/2*Sin[y[x]];
Eq2 = y''[x] == (1 - ϵ/2)*Sin[y[x]] -ϵ/2*Sin[z[x]];
bcs = {y[-xmax] == 0, y[xmax] == 2*\[Pi], z[-xmax] == 0,z[xmax] == 0};
sys = Join[{Eq1, Eq2}, bcs];
ϵ= 0.1;

NDSolve without any specified method,

s = NDSolve[sys, {y, z}, {x, -xmax, xmax}, MaxSteps -> 10^8,AccuracyGoal -> Automatic, 
PrecisionGoal -> 50,WorkingPrecision -> MachinePrecision];

NDSolve with Shooting method for StartingInitialConditions on -xmax,

soln1 = NDSolve[sys, {y, z}, {x, -xmax, xmax}, 
   AccuracyGoal -> Automatic, WorkingPrecision -> MachinePrecision, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {y'[-xmax] == 0, 
       z'[-xmax] == 0}}];

StartingInitialConditions on xmax,

soln2 = NDSolve[sys, {y, z}, {x, -xmax, xmax}, 
   AccuracyGoal -> Automatic, WorkingPrecision -> MachinePrecision, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {y'[xmax] == 0, z'[xmax] == 0}}];

StartingInitialConditions on 0,

soln3 = NDSolve[sys, {y, z}, {x, -xmax, xmax}, 
   AccuracyGoal -> Automatic, WorkingPrecision -> MachinePrecision, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {y[0] == 1.5, z[0] == 0}}];

Combining the results in one plot,

Plot[{0.5*(y[x] + z[x]) /. s, 0.5*(y[x] + z[x]) /. soln1, 
  0.5*(y[x] + z[x]) /. soln2, 0.5*(y[x] + z[x]) /. soln3}, {x, -xmax, 
  xmax}, PlotRange -> All,PlotLegends -> {"s", "soln1", "soln2", "soln3"}]

Now trying for different end points of the interval for x.

xmax=5

enter image description here

Checking whether the boundary conditions are satisfied or not?

enter image description here

xmax = 8;

enter image description here enter image description here

xmax = 15;

enter image description here

enter image description here

$\endgroup$
3
$\begingroup$

Somewhat surprisingly, I managed to find a solution with finite difference method (FDM). I've used pdetoae for discretization:

xmax = 15; ϵ = 0.1;
set = {z''[x] == (1 - ϵ/2) Sin[z[x]] - 1/2 ϵ Sin[y[x]], 
       y''[x] == (1 - ϵ/2) Sin[y[x]] - 1/2 ϵ Sin[z[x]], 
       y[-xmax] == 0, y[xmax] == 2 π, z[-xmax] == 0, z[xmax] == 0};

points = 100;
difforder = 2;
grid = Array[# &, points, {-xmax, xmax}];
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoa = pdetoae[{y, z}[x], grid, difforder];
var = # /@ grid & /@ {y, z} // Flatten;
ae = MapAt[#[[2 ;; -2]] &, ptoa@set, {{1}, {2}}];
sollst = FindRoot[ae, {#, 2 Pi} & /@ var, MaxIterations -> Infinity]; // AbsoluteTiming

Remark

If you're confused about the usage of #[[2 ;; -2]] &, the following is an alternative method for calculating sollst:

fullae = Flatten@ptoa@set;

lSSolve[obj_List, constr___, x_, opt : OptionsPattern[FindMinimum]] := 
 FindMinimum[{1/2 obj^2 // Total, constr}, x, opt]
lSSolve[obj_, rest__] := lSSolve[{obj}, rest]

sollst = lSSolve[Subtract @@@ fullae, {#, 2 Pi} & /@ var, 
     MaxIterations -> Infinity][[2]]; // AbsoluteTiming

Notice this method is slower and less robust.

sol = Partition[sollst[[All, -1]], points];
ListLinePlot[sol, PlotRange -> All, DataRange -> {-xmax, xmax}]

Mathematica graphics

The initial guess 2 Pi is important.

FindRoot still spits out lstol warning, but the solution meets the equation and boundary condition well:

{yfunc, zfunc} = ListInterpolation[#, grid] & /@ sol;
Subtract @@@ set /. {y -> yfunc, z -> zfunc};
%[[3 ;;]]
(* {0., 0., 0., 0.} *)
Plot[%%[[;; 2]] // Evaluate, {x, -xmax, xmax}, PlotRange -> All]

Mathematica graphics

The error will be smaller if you enlarge points.

I'm not sure if this is the only solution of this boundary value problem (BVP).

I believe there're better approach(es) for attacking this BVP.

$\endgroup$
2
  • $\begingroup$ Thanks dear for a nice answer. I was thinking about FDM but wasn't able to do it. $\endgroup$
    – zhk
    Jan 26, 2017 at 13:00
  • 1
    $\begingroup$ @mmm This solution is a windfall, I should say. In my (limited) experience if "Shooting" fails to solve a BVP, FDM will probably fail, too. $\endgroup$
    – xzczd
    Jan 26, 2017 at 13:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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