3
$\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?

Improve this question
$\endgroup$
4
  • $\begingroup$ Try Solve[Eq1, f''''[x]] and think about what happens at the initial condition. $\endgroup$
    – Michael E2
    May 21, 2015 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, 2015 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, 2015 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, 2015 at 15:46

2 Answers 2

Reset to default
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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Methinks one could even be bold and do f[0] == $MachineEpsilon. $\endgroup$
    – J. M.'s eventual burnout ♦
    May 21, 2015 at 20:25
2
$\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.

Improve this answer
$\endgroup$
0

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.