# NDSolve Convergence Issue for Coupled Field Problems

I have tried to use NDSolve to solve a coupled field Eqs (see the attached Codes). The MMA solver seems to failed to converge to the requested accuracy or precision within 100 iterations

ClearAll[y, x, t];
L = 10;
ode1 = y''[t] == 0;
ic11 = y == 0;
ic12 = x[L]^2*y'[L] == 100;
ode2 = -1 + x[t] + x[t]*(y'[t]^2 + y''[t]^2) - 2 x''[t] == 0;
ic21 = x' == 0;
ic22 = x'[L] == 0;
NDSolve[{ode1, ode2, ic11, ic12, ic21, ic22}, {x, y}, {t, 0, 1}]


How can I remove the error warning in MMA?

Symbolic Solution

This ODE system can be solved symbolically as follows.

sy = (DSolve[{ode1, ic11}, y, t] // Flatten) /. C -> c
(* {y -> Function[{t}, t c]} *)
sx = DSolve[{ode2 /. sy, ic21, ic22}, x, t] // Flatten
(* {x -> Function[{t}, 1/(1 + c^2)]} *)
sc = Solve[ic12 /. sx /. sy, c] // N // Flatten
(* {c -> -0.0353443 - 1.03537 I, c -> -0.0353443 + 1.03537 I,
c -> 0.0353443 - 0.964633 I, c -> 0.0353443 + 0.964633 I} *)
1/(1 + #^2) & /@ (sc // Values)
(* {-6.82769 - 7.06475 I, -6.82769 + 7.06475 I,
7.32769 + 7.06412 I,  7.32769 - 7.06412 I} *)


The fact that the solution is complex is the source of the FindRoot error.

Numerical Solution

If a numerical solution is desired, perhaps as a prototype for a more complicated system of ODEs, the following code can be used.

sn = NDSolveValue[{ode1, ode2, ic11, ic12, ic21, ic22}, {x[t], y[t]}, {t, 0, 1},
Method -> {"Shooting", "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" -> {x == -7 - 7 I, y' == -I}}];


The "Shooting" Method is needed, because this is a boundary value problem, and the "ImplicitSolver" option is needed, because the solution is complex. (The latter is illustrated here.) Note that the "StartingInitialConditions" guess does not need to be very accurate but it does need to be complex. Here are plots of the solution.

Plot[Evaluate@ReIm@First[sn], {t, 0, 1}, ImageSize -> Large,
AxesLabel -> {t, x}, LabelStyle -> {15, Bold, Black}] Plot[Evaluate@ReIm@Last[sn], {t, 0, 1}, ImageSize -> Large,
AxesLabel -> {t, y}, LabelStyle -> {15, Bold, Black}] Numerical solutions corresponding to the other values of c, above, are obtained from other choices of "StartingInitialConditions":

"StartingInitialConditions" -> {x == 7 - 7 I, y' == I}
"StartingInitialConditions" -> {x == 7 + 7 I, y' == -I}
"StartingInitialConditions" -> {x == -7 + 7 I, y' == I}


The solution above, although accurate, is incomplete in that DSolve as used omitted oscillatory eigenfunction-like solutions. They can be derived as follows:

ode2x = ode2 /. sy /. c^2 -> csq
Collect[DSolveValue[{% /. sy, ic21}, x[t], t, Assumptions -> csq < -1],
C, FullSimplify] // Flatten
(* -1 + x[t] + csq x[t] - 2 (x''[t] == 0 *)
(* 1/(1 + csq) + 2 C Cos[(Sqrt[-1 - csq] t)/Sqrt] *)


Visibly, ic22 is satisfied for n Pi/L == Sqrt[-1 - csq]/Sqrt, providing an expression forc^2 and in turn simplifying x[t].

scsq = Solve[n Pi/L == Sqrt[-1 - csq]/Sqrt, csq] // Flatten
(* {csq -> 1/50 (-50 - n^2 π^2)} *)
sn = Simplify[%% /. scsq, n > 0]
(* -(50/(n^2 π^2)) + 2 C Cos[(n π t)/10] *)


Finally, apply ic12 to evaluate C

ic12x = ic12 /. sy
Simplify[ic12x /. x -> (sn /. t -> L), n ∈ Integers];
Simplify[((#^2 & /@ %) /. c^2 -> csq) /. scsq /. C -> coef] /.
c^2 -> csq /. scsq
(* c x^2 == 100 *)
(* 1/50 (-50 - n^2 π^2) (50/(n^2 π^2) - 2 (-1)^n C)^4 == 10000 *)


From this last equation, C and in turn the final expression for x[t] are obtained, although the results are a bit long to reproduce here.

sc1 = (Solve[% /. C -> coef, coef] // Flatten) /. coef -> C
scx = sn /. # & /@ sc1


A sample plot, for two of the four n = 3 solutions, is

ReImPlot[Evaluate[scx[[3 ;; 4]] /. n -> 3], {t, 0, 10}, ImageSize -> Large, AxesLabel ->
{t, x}, LabelStyle -> {15, Bold, Black}, ReImStyle -> {Automatic, Dashed}] A corresponding NDSolve solution is

sn = NDSolveValue[{ode1, ode2, ic11, ic12, ic21, ic22}, {x[t], y[t]}, {t, 0, L},
Method -> {"Shooting", "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" -> {x == -6.5 - 5.5 I, y' == -5/3 I}}];
ReImPlot[First[sn], {t, 0, L}, ImageSize -> Large, AxesLabel -> {t, x},
LabelStyle -> {15, Bold, Black}, ReImStyle -> {Automatic, Dashed}] • thanks for your answer, Can we set "Finite Difference Method" for NDSolve? e.g. for 4th order 1D problems? Aug 11, 2020 at 11:37
• @ABCDEMMM Finite Difference is the default method for most ODEs. Finite Element also can be used in some situations, but probably not for the ODE system here. Aug 11, 2020 at 11:52
• if FD is the default method in MMA, why does NDsolve not easy solve 4th order PDEs? Aug 11, 2020 at 12:20
• @ABCDEMMM My comment above is for ODEs. PDEs are a different issue. Aug 11, 2020 at 12:33
• can NDsolve solve 4th (6th) order ODE systems? Aug 12, 2020 at 2:28