2
$\begingroup$

While solving the following ode,

Eq1 = f'''[x] - f'[x]*f'[x] + f[x]*f''[x] - K1*(2  f'[x]* f'''[x] - f[x]* f''''[x] - f''[x]*f''[x]) - f'[x] == 0

with the conditions

deqsys = {Eq1, f[0] == 0, f'[0] == 1, f'[N1] == 0, f''[N1] == 0}

and

K1 = 1.0;

Here N1=2, which is an approximation of N1=Infinity. It does not sound right but even for N1=2, I am unable to obtain a solution using NDSolve.

Now calling upon the numerical solver

NDSolve[deqsys, f, {x, 0, N1}]

I'm facing this persistent error,

Power::infy: "Infinite expression 1/0. encountered."

I have no idea, whether there is something wrong with the syntax, or the ode is a stiff one, or there is no solution to it?

$\endgroup$
  • $\begingroup$ Try Solve[Eq1, f''''[x]] and think about what happens at the initial condition. $\endgroup$ – Michael E2 May 21 '15 at 14:54
  • $\begingroup$ @Michael E2 I see your point. But when I try to solve the same equation with maple, I get solution to ode for negative values of 'K1' only. $\endgroup$ – zhk May 21 '15 at 14:59
  • $\begingroup$ Your initial conditions always contain f[0] == 0? This sort of singularity has come up on this site several times, and it each case the suggested solution is to change it to something like f[10^-6] == 10^-6. In some cases I think it should be possible mathematically to specify an initial value for the highest order derivative, but I'm not sure if NDSolve will let you. Any clue as to what step Maple takes initially? $\endgroup$ – Michael E2 May 21 '15 at 15:05
  • $\begingroup$ Well, simply using maple dsolve just handle it. As for as the background stuff is concern behind it, I don't know. $\endgroup$ – zhk May 21 '15 at 15:46
2
$\begingroup$

Here's a way that gets fairly close. Perhaps further tweaking could improve the result.

Block[{K1 = 1},
 {sol} = NDSolve[deqsys, f, {x, 0, 2}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" ->
       {f[0] == 10^-4, f'[0] == 1, f''[0] == 0.1, f'''[0] == -0.1}, 
     Method -> "StiffnessSwitching"}]
 ]

This avoids the singularity when f[x] == 0 by starting the boundary condition f[0] close to zero. Then the shooting method tries to get the initial conditions so that the boundary conditions are satisfied. They get close, but not close enough for the default NDSolve measure:

NDSolve::berr: There are significant errors {0.00267854,-4.06858*10^-11,0.0000859285,0.000110968} in the boundary value residuals. Returning the best solution found. >>

Plot[{f[x], f'[x], f''[x]} /. First[sol] // Evaluate, {x, 0, 2}]

Mathematica graphics

One can try raising "MaxIterations" higher than 100.

$\endgroup$
  • $\begingroup$ Methinks one could even be bold and do f[0] == $MachineEpsilon. $\endgroup$ – J. M. will be back soon May 21 '15 at 20:25
1
$\begingroup$

Once again, surprisingly, finite difference method (FDM) seems to help in this case. (This is against my impression that FDM will also fail if "Shooting" method fails to solve a boundary value problem. ) I'll use pdetoae for the generation of difference equation:

K1 = 1;
lb = 0; rb = 10/100;
Eq1 = f'''[x] - f'[x] f'[x] + f[x] f''[x] - 
    K1 (2 f'[x] f'''[x] - f[x] f''''[x] - f''[x] f''[x]) - f'[x] == 0;
set = {Eq1, f[lb] == 0, f'[lb] == 1, f'[rb] == 0, f''[rb] == 0};
points = 300;
difforder = 6;
grid = Array[# &, points, {lb, rb}];
(* Definition of pdetoae isn't included in this code piece,
   please find it in the link above. *)
ptoa = pdetoae[f[x], grid, difforder];
ae = MapAt[#[[3 ;; -3]] &, ptoa@set, {{1}}];
var = # /@ grid &@f // Flatten;
initialguess = 10^-2;
sollst = FindRoot[ae, {#, initialguess} & /@ var, 
    MaxIterations -> Infinity]; // AbsoluteTiming
(* {12.972383, Null} *)
func = ListInterpolation[sollst[[All, -1]], grid, InterpolationOrder -> difforder];
Plot[func[x], {x, lb, rb}, PlotRange -> All]

Mathematica graphics

Values of rb i.e. right boundary, initial guess of function value i.e. initialguess are determined by trial and error.

Error check:

Subtract @@@ set[[2 ;;]] /. f -> func
(* {0., -5.55112*10^-16, -3.017*10^-14, -7.42259*10^-10} *)
(* Residual error: *)
Plot[set[[1, 1]] /. f -> func // Evaluate, {x, lb, rb}, PlotRange -> All]

Mathematica graphics

You may feel the residual error is large, but if you enlarge points, you'll see the graph of func[x] is stable and the residual error significantly decreases.

$\endgroup$

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.