Numerical Solution of Difficult ODE Boundary Value Problems

Description

 • This page describes some strategies and suggestions for the  use of the dsolve/numeric bvp solver for difficult problems. It suggests possible solutions to be used in the case that default settings fail.
 • Initial Solution Problems
 Related errors include:
 initial Newton iteration is not converging
 initial approximate solution not converging, may need a better approximation or more mesh points
 In the case that an initial solution profile was not provided, the first error message indicates that dsolve/numeric cannot find a suitable initial solution profile. The second message indicates that because the provided solution profile is too different from the true solution, the Newton iteration is not converging.
 There are three possible strategies for this type of problem.
 1. Specify an initial approximation to the problem by using approxsoln (for the first case).
 2. Determine a "nearby" problem, and use continuation and mincont to obtain an accurate initial profile.
 3. Attempt an initial solution with a greater number of points than the default by using initmesh.
 • Continuation Problems
 This includes the errors:
 unable to refine continuation solution for parameter value ...
 continuation failed for parameter value ...
 The error messages suggest the use of a different continuation or a finer mesh, but a greater number of initial points (by using the initmesh argument) or a greater value for maxmesh may help.
 • Other Solution Errors
 Unable to Achieve Requested Accuracy
 The unable to achieve requested accuracy error indicates that either the abserr bound is set too tightly, or that maxmesh needs to be increased.
 The unable to achieve continuous solution with requested accuracy error indicates that a discrete solution with the desired error could be found, but the error estimate of the interpolant used to construct a continuous solution could not be decreased below the requested error bound. This can be fixed in the same manner as the prior error. Alternatively, it can be fixed by requesting that the interpolation error step be skipped (by specifying interpolant=false or output=mesh).
 Insufficient Precision
 The unable to achieve requested accuracy due to round-off error and precision is insufficient for required absolute error messages indicate that the resulting discretized problem is probably somewhat ill-conditioned, or that the requested tolerance is set too low. The abserr bound or Digits environmental variable must be increased. Note: The minimum number of Digits with which the solver computes is double precision (use $\mathrm{evalhf}\left(\mathrm{Digits}\right)$ to determine this number).
 Newton iteration is not converging
 This is a difficult problem. Often this error arises when a problem has narrow boundary layers or sharp corner layers. The best strategy for solving this problem is to try to obtain a low accuracy solution or a solution to a nearby problem. This may increase your understanding of the problem. Then use a change of variables to transform the boundary value problem into a more tractable problem that can be numerically solved more easily.
 • Other Errors
 Other errors can occur, such as singularity encountered, and newton iteration has not sufficiently converged after max iteration count. The first of these errors is straightforward, the solution likely has a singularity in its domain, while the second suggests that a Richardson method be used instead of a deferred correction method.
 • Singular or Stiff Problems
 The implemented methods for BVP solution are intended for non-stiff or mildly stiff BVP problems, with solutions having continuous higher order derivatives. If there are discontinuities in the higher order derivatives of the solution, the methods may be able to find a solution, but will not do so efficiently.
 For an itemized description of the options mentioned here, consult the dsolve[numeric,bvp] help page.

Examples

Example of a precision error:

 > $\mathrm{dsys}≔\left\{\mathrm{diff}\left(y\left(x\right),x,x\right)-2y\left(x\right)=0,y\left(0\right)=1.2,y\left(1\right)=0.9\right\}$
 ${\mathrm{dsys}}{≔}\left\{\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{x}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right){-}{2}{}{y}{}\left({x}\right){=}{0}{,}{y}{}\left({0}\right){=}{1.2}{,}{y}{}\left({1}\right){=}{0.9}\right\}$ (1)
 > $\mathrm{dsolve}\left(\mathrm{dsys},\mathrm{numeric},\mathrm{abserr}=1.×{10}^{-16}\right)$

Solution:

 > $\mathrm{Digits}≔20$
 ${\mathrm{Digits}}{≔}{20}$ (2)
 > $\mathrm{dsol1}≔\mathrm{dsolve}\left(\mathrm{dsys},\mathrm{numeric},\mathrm{abserr}=1.×{10}^{-16}\right)$
 ${\mathrm{dsol1}}{≔}{\mathbf{proc}}\left({\mathrm{x_bvp}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (3)
 > $\mathrm{dsol1}\left(0\right)$
 $\left[{x}{=}{0.}{,}{y}{}\left({x}\right){=}{1.2000000000000000000}{,}\frac{{ⅆ}}{{ⅆ}{x}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right){=}{-1.2525189510786390887}\right]$ (4)
 > $\mathrm{dsol1}\left(0.2\right)$
 $\left[{x}{=}{0.2}{,}{y}{}\left({x}\right){=}{0.99446362812308262645}{,}\frac{{ⅆ}}{{ⅆ}{x}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right){=}{-0.81652895731295644453}\right]$ (5)
 > $\mathrm{dsol1}\left(0.5\right)$
 $\left[{x}{=}{0.5}{,}{y}{}\left({x}\right){=}{0.83294209083370625916}{,}\frac{{ⅆ}}{{ⅆ}{x}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right){=}{-0.27638519529089584694}\right]$ (6)
 > $\mathrm{dsol1}\left(1\right)$
 $\left[{x}{=}{1.}{,}{y}{}\left({x}\right){=}{0.90000000000000000000}{,}\frac{{ⅆ}}{{ⅆ}{x}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right){=}{0.55570110924051096309}\right]$ (7)
 > $\mathrm{Digits}≔10:$

Use of Continuation:

Attempt to solve the following problem.

 > $\mathrm{ode}≔\frac{1}{100}\mathrm{diff}\left(y\left(x\right),x,x\right)+\mathrm{exp}\left(y\left(x\right)\right)\mathrm{diff}\left(y\left(x\right),x\right)-\frac{1}{2}\mathrm{\pi }\mathrm{sin}\left(\frac{1}{2}\mathrm{\pi }x\right)\mathrm{exp}\left(2y\left(x\right)\right)=0$
 ${\mathrm{ode}}{≔}\frac{\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{x}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right)}{{100}}{+}{{ⅇ}}^{{y}{}\left({x}\right)}{}\left(\frac{{ⅆ}}{{ⅆ}{x}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right)\right){-}\frac{{\mathrm{\pi }}{}{\mathrm{sin}}{}\left(\frac{{\mathrm{\pi }}{}{x}}{{2}}\right){}{{ⅇ}}^{{2}{}{y}{}\left({x}\right)}}{{2}}{=}{0}$ (8)
 > $\mathrm{bcs}≔y\left(0\right)=0,y\left(1\right)=0$
 ${\mathrm{bcs}}{≔}{y}{}\left({0}\right){=}{0}{,}{y}{}\left({1}\right){=}{0}$ (9)

An error arises indicating that the initial Newton iteration does not converge.

 > $\mathrm{dsolve}\left(\left\{\mathrm{bcs},\mathrm{ode}\right\},\mathrm{numeric}\right)$

To solve this problem, choose a continuation that increases the coefficient of the second order derivative of the equation for lambda=0, and gives the solution for lambda=1.

 > $\mathrm{newode}≔\left(\frac{1}{10}\left(1-\mathrm{\lambda }\right)+\frac{1}{100}\right)\mathrm{diff}\left(y\left(x\right),x,x\right)+\mathrm{exp}\left(y\left(x\right)\right)\mathrm{diff}\left(y\left(x\right),x\right)-\frac{1}{2}\mathrm{\pi }\mathrm{sin}\left(\frac{1}{2}\mathrm{\pi }x\right)\mathrm{exp}\left(2y\left(x\right)\right)=0$
 ${\mathrm{newode}}{≔}\left(\frac{{11}}{{100}}{-}\frac{{\mathrm{\lambda }}}{{10}}\right){}\left(\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{x}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right)\right){+}{{ⅇ}}^{{y}{}\left({x}\right)}{}\left(\frac{{ⅆ}}{{ⅆ}{x}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({x}\right)\right){-}\frac{{\mathrm{\pi }}{}{\mathrm{sin}}{}\left(\frac{{\mathrm{\pi }}{}{x}}{{2}}\right){}{{ⅇ}}^{{2}{}{y}{}\left({x}\right)}}{{2}}{=}{0}$ (10)

The 1/10 is present to make the initial solution close to the desired solution. Now apply dsolve with continuation.

 > $\mathrm{ds}≔\mathrm{dsolve}\left(\left\{\mathrm{bcs},\mathrm{newode}\right\},\mathrm{numeric},\mathrm{continuation}=\mathrm{\lambda }\right)$
 ${\mathrm{ds}}{≔}{\mathbf{proc}}\left({\mathrm{x_bvp}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (11)

The desired solution is obtained.

Use plots[odeplot](ds) to view the solution. It shows a boundary layer at x=0.

For the following problem

 > $\mathrm{ode}≔\mathrm{diff}\left(f\left(t\right),t,t,t\right)+0.5f\left(t\right)\mathrm{diff}\left(f\left(t\right),t,t\right)=0$
 ${\mathrm{ode}}{≔}\frac{{{ⅆ}}^{{3}}}{{ⅆ}{{t}}^{{3}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{f}{}\left({t}\right){+}{0.5}{}{f}{}\left({t}\right){}\left(\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{f}{}\left({t}\right)\right){=}{0}$ (12)
 > $\mathrm{bcs}≔f\left(0\right)=0,\mathrm{D}\left(f\right)\left(0\right)=0,\mathrm{D}\left(f\right)\left(100\right)=1$
 ${\mathrm{bcs}}{≔}{f}{}\left({0}\right){=}{0}{,}{\mathrm{D}}{}\left({f}\right){}\left({0}\right){=}{0}{,}{\mathrm{D}}{}\left({f}\right){}\left({100}\right){=}{1}$ (13)

you cannot obtain the solution directly:

 > $\mathrm{dsolve}\left(\left\{\mathrm{bcs},\mathrm{ode}\right\},\mathrm{numeric}\right)$

There is an easy solution of $f\left(x\right)=0$ if the right boundary condition was $\mathrm{D}\left(f\right)\left(100\right)=0$, so use this as the continuation, and obtain a solution:

 > $\mathrm{newbcs}≔f\left(0\right)=0,\mathrm{D}\left(f\right)\left(0\right)=0,\mathrm{D}\left(f\right)\left(100\right)=\mathrm{\lambda }$
 ${\mathrm{newbcs}}{≔}{f}{}\left({0}\right){=}{0}{,}{\mathrm{D}}{}\left({f}\right){}\left({0}\right){=}{0}{,}{\mathrm{D}}{}\left({f}\right){}\left({100}\right){=}{\mathrm{\lambda }}$ (14)
 > $\mathrm{dsn}≔\mathrm{dsolve}\left(\left\{\mathrm{newbcs},\mathrm{ode}\right\},\mathrm{numeric},\mathrm{continuation}=\mathrm{\lambda }\right)$
 ${\mathrm{dsn}}{≔}{\mathbf{proc}}\left({\mathrm{x_bvp}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (15)

and with the following plot, there is a boundary layer in the second derivative of the solution at t=0

 > $\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsn},\left[t,\mathrm{diff}\left(f\left(t\right),t,t\right)\right]\right)$