# Using shooting method to solve coupled ODE for domain walls

How can one solve a coupled equation for f[x] and a[x] with known boundary conditions for f[x] and a'[x] at some -xmax and xmax?

κ=0.2; xmax=10;

eq1={1/κ^2 f''[x] == (a[x]^2 - 1) f[x] + f[x]^3};
eq2={a''[x] == a[x] f[x]^2};
boundaryconditions={f[-xmax] == 0, f[xmax] == 1, a'[-xmax] == 1/Sqrt, a'[xmax] == 0};


An attempt to use NDSolve with the shooting method didn't work.

## 1 Answer

For this particular problem, NDSolve with Method -> {"Shooting", "StartingInitialConditions" -> ...}} becomes ever more sensitive to the StartingInitialConditions selected as xmax is increased. It is, in essence, a separatrix problem. I used a brute force approach to reach xmax == 76/10:

κ = 1/5; xmax = 76/10;
eq1 = {1/κ^2 f''[x] == (a[x]^2 - 1) f[x] + f[x]^3};
eq2 = {a''[x] == a[x] f[x]^2};
bc = {f[-xmax] == 0, f[xmax] == 1, a'[-xmax] == 1/Sqrt, a'[xmax] == 0};
sol = NDSolveValue[{eq1, eq2, bc}, {f[x], a[x]}, {x, -xmax, xmax},
Method -> {"Shooting", "StartingInitialConditions" -> {f'[-xmax] == 2883 10^-5,
a[-xmax] == -4324 10^-3}}, WorkingPrecision -> 30, PrecisionGoal -> 10] // Flatten;
Plot[sol, {x, -xmax, xmax}, AxesLabel -> {x, "f, a"},
LabelStyle -> Directive[Black, Bold, Medium], ImageSize -> Large] {D[sol[], x], sol[]} /. x -> -xmax
(* {0.028833219726590837125, -4.3237695236626951946629746429} *)


The approach used was to solve the problem for xmax = 2, for which the NDSolve automatic shooting works well. Then, determine {D[sol[], x], sol[]} /. x -> -xmax from that solution to obtain a StartingInitialConditions guess for a larger value of xmax, and so on. At first, the incremental increase in xmax can be fairly large, but eventually it must be less than one part in 100.

Improved Result

As discussed in comments below, Shooting from the center of the domain (here, x == 0) cuts the maximum integration distance in half, although at the cost of requiring FindRoot, called internally by NDSolve, to solve for four constants instead of two. For the present problem, Shooting from x == 0 works well in that the iterative approach described above requires many fewer steps.

κ = 1/5; xmax = 10;
eq1 = {1/κ^2 f''[x] == (a[x]^2 - 1) f[x] + f[x]^3};
eq2 = {a''[x] == a[x] f[x]^2};
bc = {f[-xmax] == 0, f[xmax] == 1, a'[-xmax] == 1/Sqrt, a'[xmax] == 0};
soc = NDSolveValue[{eq1, eq2, bc}, {f[x], a[x]}, {x, -xmax, xmax},
Method -> {"Shooting", "StartingInitialConditions" -> {f == 7060 10^-4,
a == -1048 10^-4, f' == 0716 10^-4, a' == 0784 10^-4}},
WorkingPrecision -> 30, PrecisionGoal -> 10] // Flatten;
Plot[soc, {x, -xmax, xmax}, AxesLabel -> {x, "f, a"}, PlotRange -> All,
LabelStyle -> Directive[Black, Bold, Medium], ImageSize -> Large] {soc, D[soc, x]} /. x -> 0
(* {{0.709250089398753089314125661908, -0.101241412225167625937997398726},
{0.0709665819441012826702016091966, 0.0759572604513366536795123139364}} *)

• The high level of sensitivity to starting initial conditions seems to be a common issue when working with equations of motion relevant to Ginzburg-Landau theory. The issue is that the shooting method has to shoot to a specified boundary condition, while also making the derivative at the boundary close to zero with some tolerance. – 121 May 2 '18 at 18:20
• An approach that may be beneficial is to shoot from both ends with a very shallow gradient close to the boundaries and then match the solutions in the middle with a specified tolerance. – 121 May 2 '18 at 18:27
• Shallow gradient for both f[x] and a'[x], which is the magnetic field in the domain wall problem. – 121 May 2 '18 at 18:36
• @121 NDSolve can integrate from any point in the domain to both endpoints. Shooting from the center has the advantage of cutting the maximum integration distance in half but has the disadvantage of requiring FindRoot, which is called by NDSolve when using Shooting, to solve for four constants instead of two. – bbgodfrey May 2 '18 at 18:37
• @bbgodfrey I have a similar but different question. Could you take a look? Thanks! – Boson Bear May 20 '18 at 3:22