Solution of a two point boundary problem: finite differences

directly to the input form: predefined cases

directly to the input form: user defined case

• In the following, y, f and g₁, g₂ denote real functions with one, three respectively two variables. The prime denotes differentiation with respect to a ''free'' real variable x. This ''free'' real parameter x has usually here the physical meanig of ''space'' (and only seldom ''time''). We search for a solution y as function of x, given some relation between y, y', y'' and x and some conditions on y, y' at two points x=a, x=b, the so called boundary conditions.
  Specifically the problem dealt with here is:
  Find an (at least) two times differentiable function y defined on the interval [a,b] such that there holds
  y''(x) = f(x,y(x),y'(x)) for a < x < b
  and
  g₁(y(a),y'(a)) = 0 , g₂(y(b),y'(b)) = 0 .
  That means that we restrict the problem to so called separated boundary conditions. The theory also considers the case where g₁ and g₂ both depend on the values taken at a and b.
• The essential difference opposed to initial value problems is the fact that the length of the interval of existence of the solution is part of the problem setting and that conditions on the solution are given at two different points. This type of problem is called a boundary value problem, BVP for short. For the case of f and g₁, g₂ linear in y, y', there exists a complete existence and uniqueness theory, but for the general nonlinear case no such general results exist.
• We present here the simplest numerical approach: the so called finite difference method which consists in three steps:
  1. One considers the differential equation on a finite grid in [a,b] only.
  2. One replaces the values of the derivatives of y at the grid points by derivatives of interpolation polynomials interpolating the values of y at neighboring grid points.
  3. This results in a finite system of (non)linear equations for the unknown values of y at the grid points which is then solved by some standard solver like those presented here in the chapter on linear and nonlinear equations.
• The finite differences used here are the simplest ones and we restrict us here to an equidistant grid with grid size h > 0.
  1. (y_i+1-y_i)/h for y'_i or y'_i+1 at the boundaries and
     (y_i+1-y_i-1)/(2h) for y'_i in the interior of [a,b]
  2. (y_i+1-2y_i+y_i-1)/h² for y''_i
• This results in global discretization errors which tend to zero like Ch² if there are only Dirichlet (function values) boundary conditions and otherwise the convergence is only like Ch, a precision much less than we are used to employ in solving initial value problems. In many cases however we may rise the order of convergence by the technique of Richardson extrapolation in the same manner like we use it in numerical differentiation and in Rombergs method of quadrature. If we can assume an expansion in even powers of h only, we use order 2 extrapolation leading to order 6 and in the other case we use an extrapolation scheme covering also odd orders of h leading to order 3 (at most). For some more details see here .
• In order to use a better approximation to the derivative in the boundary conditions there is a technique known as the use of ''fictitious points'', which we provide too. One can use it if the functions f, g₁, g₂ are well defined on an interval containing [a,b] in its interior. In this case the convergence is always of second order in h.
• We always try to solve the problem on three coherent grids with grid sizes h, h/2, h/4 in order to let you check the convergence order and the effect of Richardson extrapolation.

There are two different input forms: Predefined
and User defined

Input: Predefined cases

• There is a choice between 6 different problems, of which the analytical solution is known. You find the complete description of these cases here . You select a problem by clicking on the corresponding box.
• In case of problem 3, which has two solutions, you can select which one is desired.
• Next you decide whether to use the technique of fictitious points or not.
• You select the number of grid points N+1 by selecting N. The grid size is then
    h = (b-a)/N
  
• Next you decide whether to use the predefined internal parameters
  itmax, the maximum number of Newton steps in the solution process of the discrete (non)linear system of equations, the required relative precision eps for the solution of this system (with the absolute precision set to eps²) and printtable which gives you the possibility to get a printed table of the solution on the crudest grid, or to set these values as desired. The default values are 30, 10^-10, false .
• You have the possibility to set the initial guess for the solution of the system of equations yourself or to generate this internally. In the latter case this initial guess is generated from the known analytical solution by disturbing it by pseudo random errors controlled by a variable perturb. The largeness of the perturbations is perturb × max(abs(y(x_i))). In the case of a linear problem (case 1) this plays no role since a linear system is solved by Newtons method in one step irrespective of the initial guess, but in the other nonlinear cases a large perturbation may cause considerable trouble.

Output: predefined case

• The final boundary values for the three grids and the corresponding extrapolated value.
• The required computing time.
• If requested a printed table of the computed solution on the coarsest grid.
• (Up to) five plots:
  1. the discretization errors on the three grids and for the extrapolated values.
  2. the development of the residual norm during the Newton process. The three curves correspond to the three grids in the sequence red - green -blue from the crudest to the finest
• Should the BVP-solver issue warning messages you get these as a protocol.

Questions: Predefined case?!

• Behaves the discretization error of the solver as expected ?
• In which cases fails the Richardson extrapolation and why?
• Which effect is to be seen if you take a too large h?

to the input form: predefined case

Input: user defined case

• Here you have the possibility to define a BVP yourself. This requires
  1. Specification whether to use fictitious points or not.
  2. Specification of a,b and of N with
       h = (b-a)/N
     
  3. Whether to see a printed table of the solution.
  4. You write a piece of code for computing f and two more pieces of code for computing g₁ and g₂. Details are given on the input form and in the help text.
  5. If you know the analytical solution of your problem, then you write a fourth piece of code for computing it. In this case you may generate the initial guess for the solution from the exact solution by disturbing it by pseudo random errors controlled by a variable perturb. The largeness of the perturbations is perturb × max(abs(y(x_i))).
  6. If you know no analytical solution or don't want to use it you must specify the N+1 resp. N+3 initial values yourself.

Output: user defined case

• This is the same as for the predefined cases with the exception that in the case of no known analytical solution the discretized solution is shown instead of the discretization error.

Questions: user defined case?!

• This case is meant for the solution of a BVP problem of your own interest.
• The questions to be asked here are of course the same as for the predefined cases.

to the input form: user defined case