# NDSolve issuing boundary condition error message

I would like to solve this non-linear differential equation with boundary conditions and a small perturbation (eps):

eps = 0.01;
eq = (eps^2) y''[x] + eps*x*y'[x] - y[x] == -Exp[x];
bcs = {y[0] == 2, y[1] == 1};
s = NDSolve[{eq, bcs}, y, {x, 0, 1}]
Plot[Evaluate[y[x] /. s], {x, 0, 1}, PlotRange -> All]


But the solution does not fit the BC: (expected due to this error)

NDSolve: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly.

So I tried different methods such as : "Chasing" or "Shooting". (Infos), but none of them match the both bondary conditions simultaneously.

Any ideas?

• It seems like you're perturbing the equation y[x] == Exp[x], which itself does not satisfy the BCs. It's not surprising it's ill-conditioned. It's probably unstable as eps goes to zero. I suppose this is what you're studying? Mar 17, 2018 at 16:15

## Introduction

I would have done a Chebyshev approach earlier, but a stupid, elusive typo kept the accuracy low. The method is fast and highly accurate, although for the code below, recursive refinement and a convergence test have not been implemented (the idea would be the same as used for chebapprox in this answer).

A discussion of this approach, as well as an example over an infinite interval, may be found in this answer. Here we need only a simple affine transformation of a finite interval to the standard Chebyshev interval $[-1, 1]$. See Boyd, Chebyshev and Fourier Spectral Methods (2001) for more.

## OP's problem

For convenience, here is the OP's setup:

eps = 1/100;
eq = (eps^2) y''[x] + eps*x*y'[x] - y[x] == -Exp[x];
bcs = {y[0] == 2, y[1] == 1};


### Chebyshev series solutions of linear, second-order BVPs

First we give a general function that will solve a linear, second-order BVP. It depends on some utilities given at the end.

ClearAll[chebLinearSecondOrderBVP, chebLinearSecondOrderBVPparse];
Options[chebLinearSecondOrderBVP] = {WorkingPrecision -> MachinePrecision};
SetAttributes[chebLinearSecondOrderBVPparse, HoldAllComplete];
(* parses the input; TBD: error messages! *)
chebLinearSecondOrderBVPparse[
chebLinearSecondOrderBVP[sys_, y_, {x_, a_?NumericQ, b_?NumericQ},
n_, opts : OptionsPattern[chebLinearSecondOrderBVP]]] :=
Module[{ode, bcs, odebcs, dorder, y1, y2, prec},
With[{f = FreeQ[#, y[x] | Derivative[_][y][x]] &},
odebcs = SortBy[f]@GatherBy[sys, f]
];
If[Length /@ odebcs == {1, 2},
{{ode}, bcs} = odebcs,
Return[$Failed]]; {y1, y2} = {y[a], y[b]} /. First@Solve[bcs, {y[a], y[b]}]; If[! VectorQ[{y1, y2}, NumericQ], Return[$Failed]];
dorder = Max@Cases[ode, Derivative[k_][y][x] :> k, Infinity];
If[dorder != 2, Return[$Failed]]; If[! InternalLinearQ[ode, {y[x], y'[x], y''[x]}], Return[$Failed]];
prec = OptionValue[chebLinearSecondOrderBVP, opts, WorkingPrecision];
{ode, y1, y2, prec}
];
call : chebLinearSecondOrderBVP[
sys_, y_, {x_, a_?NumericQ, b_?NumericQ}, n_, OptionsPattern[]] :=
Module[{ode, y1, y2, xvec, dm, dops, opL, load, dy, yvec, cc, prec, res},
res = chebLinearSecondOrderBVPparse[call]; (* parse arguments *)
(
{ode, y1, y2, prec} = res;
(* construct differential operator  opL  for the ode *)
ode = Normal@ CoefficientArrays[ode /. Equal -> Subtract, {y[x], y'[x], y''[x]}];
{dm, {{xvec}}} = Reap[chebDM[n, prec], "x"];  (* Ch. deriv. operator and pts. *)
dops = NestList[2/(b - a) dm.# &, IdentityMatrix[n + 1], 2]; (* 0-2 order diff.ops.*)
dops = dops[[All, 2 ;; -2, 2 ;; -2]];         (* strip boundary *)
xvec = Rescale[xvec[[2 ;; -2]], {-1, 1}, {a, b}];
opL = (ode[[2]] /. x -> xvec) dops // Total;  (* assumes coefficients are Listable *)
(* construct load and subtract line to make homogeneous BCs y[a]==y[b]==0:
* the line affects the y[x], y'[x] terms (ode[[2,1;;2]]);
* it does not affect the y''[x] term
*)
load = ode[[1]] /. x -> xvec;
dy = {(a y2 - b y1 + x (y1 - y2))/(a - b), (y2 - y1)/(b - a)} ode[[2, 1 ;; 2]] /.
x -> xvec // Total;                        (* impose homogeneous BCs *)
(* solve BVP, get Chebyshev coeffs and add line back for BVs *)
yvec = LinearSolve[opL, -load - dy];
cc = chebSeries[yvec];
cc[[1]] += (y1 + y2)/2;  (* Adjust to BCs *)
cc[[2]] += (y2 - y1)/2;
res = chebFunc[cc, {a, b}]
) /; FreeQ[res, \$Failed | chebLinearSecondOrderBVPparse]
];


The solution of the OP's problem:

ySol2 = chebLinearSecondOrderBVP[{eq, y[0] == 2, y[1] == 1},
y, {x, 0, 1}, 64, WorkingPrecision -> 32];

Plot[{Exp[x], ySol2[x]}, {x, 0, 1}]


The solution satisfies the BVP to machine precision

Plot[Evaluate[eq /. Equal -> Subtract /. y -> ySol2], {x, 0, 1},
PlotRange -> All, PlotPoints -> 200]


{y[0.] - 2, y[1.] - 1} /. y -> ySol2  (* BCs *)
(*  {0., 0.}  *)


### Chebyshev utilities

(*
* Chebyshev series utilities
*)

ClearAll[chebFunc];
chebFunc::usage = "f = chebFunc[c,{a,b}], c = {c0, c1,..., cn} Chebyshev coefficients, over the interval {a,b}; y = chebFunc[c,{a,b}][x] evaluates the function";
chebFunc[c_, {a_, b_}][x_] := chebFunc[c, {a, b}, x];
chebFunc[c_?(VectorQ[#, NumericQ] &), {a_?NumericQ, b_?NumericQ}, x_?NumericQ] :=
ChebyshevT[Range[0, Length[c] - 1], (2 x - (a + b))/(b - a)].c;
chebFunc[c_?(VectorQ[#, NumericQ] &), {a_?NumericQ, b_?NumericQ}, x_?(VectorQ[#, NumericQ] &)] :=
Cos[Outer[Times, ArcCos[(2 x - (a + b))/(b - a)], Range[0, Length[c] - 1]]].c;
chebFunc /: Normal[chebFunc[c_?VectorQ, {a_, b_}, x_]] :=
Evaluate@ChebyshevT[#, (2 x - (a + b))/(b - a)] &[Range[0, Length[c] - 1]].c;
Derivative[_, _, n_][chebFunc][c_, {a_, b_}, x_] /; 0 <= n <= Length@c :=
chebFunc[Nest[chebDerivative, c, n] (2/(b - a))^n, {a, b}, x];

ClearAll[chebDerivative, chebSeries, chebDM];
chebDerivative::usage = "chebDerivative[c, {a,b}] differentiates the Chebyshev series c scaled over the interval {a,b}";
chebDerivative[c_] := chebDerivative[c, {-1, 1}];
chebDerivative[c_, {a_, b_}] :=
Module[{c1 = 0, c2 = 0, c3},
2/(b - a) MapAt[#/2 &,
Reverse@Table[
c3 = c2;
c2 = c1;
c1 = 2 (n + 1)*c[[n + 2]] + c3,
{n, Length[c] - 2, 0, -1}],
1]];
chebDM::usage = "chebD[n] = (n+1)x(n+1) Chebyshev spectral differentiation matrix of order n";
chebDM[n_Integer?Positive, prec_ : MachinePrecision] :=
Module[{fac, xm},
xm = If[prec === MachinePrecision,
ConstantArray[Sin[(π Range[N[n], -N[n], -2.])/(2 n)], n + 1],
ConstantArray[N[Sin[(π Range[n, -n, -2])/(2 n)], prec], n + 1]];
Sow[First[xm], "x"];
fac = Flatten[{2, PadRight[{}, n - 1, {-1, 1}], 2 - 4 Mod[n, 2]}];
xm = Outer[Times, fac, 1/fac]/(Transpose[xm] - xm + IdentityMatrix[n + 1]);
xm - DiagonalMatrix[Total[xm, {2}]]];

chebSeries::usage = "chebSeries[y] returns the Chebyshev series for the function values y over the standard interval {-1, 1}";
chebSeries[y_] := Module[{cc},
cc = Sqrt[2/(Length[y] - 1)] FourierDCT[y, 1];
cc[[{1, -1}]] /= 2;
cc];


Here are a couple more ways. The first is similar to @anderstood's, but with a different integration method. The second uses FEM and has the virtue of being fast. It appears to be a perturbation of the equation y[x] == Exp[x] with discontinuous boundary conditions y[0] == 2, y[1] == 1. The problem suffers some from boundary effects. The options for the first method are set so that this virtually eliminated, except slightly near x == 1. The PrecisionGoal has to be set fairly high or the ill-conditioning to great an error. But once it's over the threshold, you get a highly accurate solution. The FEM is limited to machine precision, and it effectively achieves PrecisionGoal -> 8 when far enough away from the boundary. If high precision is not needed, it's probably the best method. Note you have to set the "MaxCellMeasure" for it depending on eps, perhaps eps/10 is good enough.

Clear[x, y];
eps = 1/100;
eq = (eps^2) y''[x] + eps*x*y'[x] - y[x] == -Exp[x];
bcs = {y[0] == 2, y[1] == 1};

s = NDSolve[{eq, bcs}, y, {x, 0, 1}, Method -> "Extrapolation",
PrecisionGoal -> 32, WorkingPrecision -> 60,
InterpolationOrder -> All];

s2 = NDSolve[{eq, bcs}, y, {x, 0, 1},
Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 1/1000}}];


Both produces indistinguishable graphs. I compare one with the graph of Exp[x]:

Plot[Evaluate[{Exp[x], y[x] /. s}], {x, 0, 1}, PlotRange -> All]


Evaluating the residual of the ODE shows the boundary effects:

Plot[Evaluate[
eq /. Equal -> Subtract /. Join[s, s2] // RealExponent],
{x, 0, 1}, PlotRange -> {-36, 0}, PlotPoints -> 200,
WorkingPrecision -> 32]


Increase the working precision (see 22081) plus "StiffnessSwitching" method:

eps = 1/100;
eq = (eps^2) y''[x] + eps*x*y'[x] - y[x] == -Exp[x];
bcs = {y[0] == 2, y[1] == 1};
s = NDSolve[{eq, bcs}, y, {x, 0, 1}, Method -> "StiffnessSwitching",
WorkingPrecision -> 80]
Plot[Evaluate[y[x] /. s], {x, 0, 1}, PlotRange -> All]
`