I tried to solve this differential equation:

$$\epsilon y''(x)+xy'(x)=-\epsilon \pi^2 \cos(\pi x)-\pi x\sin(\pi x)$$

with boundary conditions: $y(-1)=-2, \space y(1)=0$. If we take $\epsilon=0.1$, Mathematica can solve it without any trouble

Block[{e = 0.1, min = -1, max = 1},
   y[x] /. NDSolve[{e y''[x] + 
        y'[x] x == -e Pi^2 Cos[Pi x] - Pi x Sin[Pi x], 
      y[min] == -2, y[max] == 0}, y, {x, min, max}]], {x, min, max}]

enter image description here

But if we want a smaller $\epsilon$, let say 0.01, Mathematica seems unable to handle it. Is there any options to invoke or methods to employ to get the desired result? Anyway, this is the solution for $\epsilon=0.0001$.

enter image description here

Thank you.

  • $\begingroup$ Piece of advice: don't use Block to inject values into parameters. Use With instead. $\endgroup$
    – m_goldberg
    Commented Jun 25, 2015 at 14:24
  • 2
    $\begingroup$ The fundamental problem would appear to be that, in the limit of small e, the order of the equation drops from second to first, with the result that there is one too many boundary conditions. $\endgroup$
    – bbgodfrey
    Commented Jun 25, 2015 at 14:30

3 Answers 3


DSolve can handle this.

y[x_, e_] = y[x] /. DSolve[{
      e y''[x] + y'[x] x == -e Pi^2 Cos[Pi x] - Pi x Sin[Pi x],
      y[-1] == -2, y[1] == 0}, y[x], x][[1]] // Simplify

enter image description here

 Plot[y[x, e], {x, -1, 1}],
 {{e, 0.01}, 0.0001, 0.1, Appearance -> "Labeled"}]

enter image description here

  • $\begingroup$ Thank you for your answer, but I'm looking for numerical solution. $\endgroup$
    – Deco
    Commented Jun 25, 2015 at 14:59
  • 2
    $\begingroup$ @Deco Bob's symbolic solution can be evaluated numerically quickly, and with better accuracy than a mere numerical approximation. In fact DSolve can solve the problem in full generality: DSolve[{e y''[x] + y'[x] x == -e Pi^2 Cos[Pi x] - Pi x Sin[Pi x], y[min] == -2, y[max] == 0}, y, x]. Is there some reason you do not think this is a superior solution? $\endgroup$
    – Michael E2
    Commented Jun 25, 2015 at 15:55

As noted by @bbgodfrey, the "shooting" algorithms that Mathematica tends to use are not well-adapted to this particular equation. Better would be some kind of relaxation method, which is what Mathematica uses (I think) for solving PDEs on a mesh. And an ODE is just a PDE in one dimension, so let's try solving this equation on a one-dimensional mesh:

truey[x_,e_] = Cos[\[Pi] x] + Erf[x/(Sqrt[2] Sqrt[e])]/Erf[1/(Sqrt[2] Sqrt[e])];

e = 0.0001; min = -1; max = 1;
mesh = ToElementMesh[Interval[{min, max}], MaxCellMeasure -> 0.05];
bcs = {DirichletCondition[y[x] == -2, x == min], DirichletCondition[y[x] == 0, x == max]}
soln = NDSolve[{e y''[x] + y'[x] x == -e Pi^2 Cos[Pi x] - Pi x Sin[Pi x], bcs}, y, Element[x, mesh]]
Plot[{Evaluate[y[x] /. soln], truey[x, e]}, {x, min, max}]
Clear[e, min, max]

enter image description here

I've plotted here the result from Mathematica's finite element solver (in blue) versus the true solution found by @BobHanlon (in yellow). I've actually used a coarser mesh than is ideal to show the difference between the two; if you set MaxCellMeasure -> 0.01 in the above code (instead of 0.05), the two curves are indistinguishable at this resolution.


To see why NDSolve has difficulty with this problem for very small e, consider that NDSolve solves this two-point boundary value problem by some form of shooting. In other words, it varies y'[min] until one is found that yields the desired y[max]. However, as e becomes very small, the sensitivity of y[max] to y'[min] becomes great, because the differential equation becomes singular in that limit.

y''[x] + y'[x] x/e == -Pi^2 Cos[Pi x] - Pi x Sin[Pi x]/e

This is also apparent from Bob Hanlon's symbolic solution. To illustrate the sensitivity of y[max] to y'[min], consider

slope = D[Cos[π x] + Erf[x/(Sqrt[2] Sqrt[e])]/Erf[1/(Sqrt[2] Sqrt[e])],  x] /. x -> -1;
LogLogPlot[N[slope], {e, 0.0001, .1}, PlotRange -> All]

enter image description here

You can work around this sensitivity by increasing WorkingPrecision, e. g. for e = .01,

With[{e = 1/100, min = -1, max = 1}, 
    Plot[Evaluate[y[x] /. NDSolve[{e y''[x] + y'[x] x == -e Pi^2 Cos[Pi x] - 
      Pi x Sin[Pi x], y[min] == -2, y[max] == 0}, y, {x, min, max}, 
      WorkingPrecision -> 50, MaxSteps -> 50000]], {x, min, max}]]

enter image description here

but doing so rapidly becomes prohibitively expensive as e is further reduced.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.