Numerical solution of boundary value problems

Linear, second order

Although the standard Mathematica routine NDSolve is supposed to work for this class of problems, it is remarkable delicate, and fails often even if ϵ is not particularly small. It uses the idea of 'shooting'. Using for ex. the left boundary condition (BC) and guessing on a second left condition (say, for the first derivative) two different solutions are found numerically. Typically neither of these two will satisfy the right BC. However, because of linearity, one can explicitly find a linear conbination of the two trial solutions will do this.
		
We use below Mathematica's built-in solver for the example in B&O page 431.

ode = 
      ϵ y''[x] + (1 + x) y'[x] + y[x] == 0

sol = 
      NDSolve[{ode  /. ϵ -> 0.1, y[0] == 1, y[1] == 1}, 
        y[x], {x, 0, 1}]

The output of NDSolve takes the form of an InterpolatingFunction - essentially a table of intermediate values from which the function values at arbitrary points can be obtained by standard interpolation. We plot the result as follows

Plot[Evaluate[y[x]  /.  sol], {x, 0, 1}]

Mathematica fails for this problem if we try ϵ=0.05 - still not a very small value.

Nonlinear, or higher order equations

We can now implement basically the same principle ourselves. Instead of just using two tries, we  have to allow several tries, and use Mathematica's root finding algorithm (for general nonlinear functions)  in order to converge in on the solution. 

As an example of this, we consider a third order ODE from B&O page 446 (linear, but third order; Mathematica will not accept the direct approach above):

ode  =  ϵ y'''[x] - y'[x] + x y[x] == 0

with the BC  y[0]=y'[0]=y[1]=1. Since two conditions are already given at x=0, we decide to there add a third; y''[0]=a. The problem has now become to choose the value of a so that y[1]=1. Given any value of a, we can compute the corresponding value of y[1] through the following function

\(r[a_] := \ 
      \((sol = NDSolve[{ode\  /. \ ϵ -> 0.01, y[0] == 1, 
                \(y'\)[0] == 1, \(\(y'\)'\)[0] == a}, y[x], {x, 0, 1}]; 		
          Evaluate[y[x] /. sol]\  /. x -> 1)\)\ [\([1]\)]\)

We could for example try a few values of the variable a 'manually' to see if we can get the result to become close to 1:

r[0]

r[-10]

r[-9]

r[-8.9]

Mathematica has a function FindRoot that does this kind of search for us. After specifying the equation to solve, FindRoot requires one or two starting guesses, and then conducts a local search. With one guess, it tries a Newton-type method - and fails in this type of cases. With two initial guesses (which normally don't have to be very good at all) it successfully tries a secant method. So

FindRoot[r[a] == 1, {a, 0, 1}]

FindRoot called r[a] several times, and we can still access the full solution generated by its last ODE evaluation. So we plot this

Plot[Evaluate[y[x]  /.  sol], {x, 0, 1}]