2
$\begingroup$

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] == 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] == 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?

$\endgroup$
6
$\begingroup$

Symbolic Solution

This ODE system can be solved symbolically as follows.

sy = (DSolve[{ode1, ic11}, y, t] // Flatten) /. C[2] -> 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[0] == -7 - 7 I, y'[0] == -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}]

enter image description here

Plot[Evaluate@ReIm@Last[sn], {t, 0, 1}, ImageSize -> Large, 
    AxesLabel -> {t, y}, LabelStyle -> {15, Bold, Black}]

enter image description here

Numerical solutions corresponding to the other values of c, above, are obtained from other choices of "StartingInitialConditions":

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

Addendum: Oscillatory Solutions

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[1], FullSimplify] // Flatten
(* -1 + x[t] + csq x[t] - 2 (x''[t] == 0 *)
(* 1/(1 + csq) + 2 C[1] Cos[(Sqrt[-1 - csq] t)/Sqrt[2]] *)

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

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

Finally, apply ic12 to evaluate C[1]

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

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

sc1 = (Solve[% /. C[1] -> coef, coef] // Flatten) /. coef -> C[1]
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}]

enter image description here

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[0] == -6.5 - 5.5 I, y'[0] == -5/3 I}}];
ReImPlot[First[sn], {t, 0, L}, ImageSize -> Large, AxesLabel -> {t, x}, 
    LabelStyle -> {15, Bold, Black}, ReImStyle -> {Automatic, Dashed}]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ thanks for your answer, Can we set "Finite Difference Method" for NDSolve? e.g. for 4th order 1D problems? $\endgroup$ – ABCDEMMM Aug 11 at 11:37
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Aug 11 at 11:52
  • $\begingroup$ if FD is the default method in MMA, why does NDsolve not easy solve 4th order PDEs? $\endgroup$ – ABCDEMMM Aug 11 at 12:20
  • $\begingroup$ @ABCDEMMM My comment above is for ODEs. PDEs are a different issue. $\endgroup$ – bbgodfrey Aug 11 at 12:33
  • $\begingroup$ can NDsolve solve 4th (6th) order ODE systems? $\endgroup$ – ABCDEMMM Aug 12 at 2:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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