Solution of a two point boundary problem: local collocation

directly to the input form: predefined cases

directly to the input form: user defined case

• In the following, y(x), F(x,y) and R(ya,yb) denote real vector functions with the real variable x and the vector variables y, ya, yb of dimension n. The prime denotes differentiation with respect to the ''free'' real variable x. This ''free'' real variable 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' and x and some conditions on y at two points x=a, x=b, the so called boundary conditions.
  Specifically the problem dealt with here is:
  Find an (at least) one times differentiable function y defined on the interval [a,b] such that there holds
  y'(x) = F(x,y(x)) for a < x < b
  and
  R(y(a),y(b)) = 0 .
  
• 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 R 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 currently mostly used technique for its solution, the method of local collocation.
  This method consists in three steps:
  1. One subdivides [a,b] in intervals [x_i,x_i+1], i=1,..,N-1 with x₁ = a , x_N = b.
  2. One replaces the function y on each subinterval by some ''simple'' function fixed by a number of parameters, e.g. a polynomial or a spline. This function will be uniquely defined by the requirement taking at x_i the value of the corresponding piece at the left (resp. ya for i=1) and satisfying the differential function at given nodes inside [x_i,x_i+1].
  3. In addition one requires that the function composed from these pieces satisfies the boundary conditions.
  For example in the simplest case used here we choose a vector of linear polynomials fixed by their value at x_i and satisfying the differential equation in the middle of the interval: p_i(x)=p_i,0 + (x - x_i)*p_i,1 , i.e. p_i'(x) = p_i,1
  
  1. p_i(x_i) = y_i , fixing p_i,0
  2. p'_i = F ( (x_i+x_i+1)/2 , p((x_i+x_i+1)/2) ), fixing p_i,1
  and set y_i+1 = p_i(x_i+1) , which makes the composed function continuous.
  
• This method is called the implicit midpoint rule.
• This results in a finite system of (non)linear equations for the unknown values of y_i at the grid points and can be rewritten as
  
  1. y_i+1 - y_i -(x_i+1 - x_i)F((x_i+x_i+1)/2,(y_i+1 + y_i)/2) = 0 , i=1,..,N-1
  2. R(y₁,y_N) = 0 .
• This system is then solved by some standard solver like those presented here in the chapter on linear and nonlinear equations. We use the damped Newton method, which evaluates the Jacobian in each step. The Jacobian of this system has a special sparse structure

  A1,1	A1,2	O	O	...	O
  O	A2,2	A2,3	O	...	O
  ...	...	...	...	...	...
  O	...	...	O	AN-1,N-1	AN-1,N
  JRya	O	...	...	O	JRyb

  Here
  A_i,i = -I - (h_i/2) JF((x_i+x_i+1)/2,(y_i+1+y_i)/2))
  A_i,i+1 = I - (h_i/2) JF((x_i+x_i+1)/2,(y_i+1+y_i)/2)) ,
  h_i = x_i+1 - x_i .
  JF denotes the Jacobian of F w.r.t. y and JR_ya,yb denotes the Jacobian of R w.r.t. to its two arguments ya, yb. We take advantage of this structure by decomposing the first N-1 block rows separately and exclude thereby the last block row from pivoting which comes into play only in decomposing the last block row. This can introduce pivotal growth which therefore is reported in the results and can be reduced by decreasing the stepsize h_i (here by using a finer initial grid for example). However, for problems with boundary singularities this might fail.
• This method converges for two times differentiable F and R (w.r.t. to x,y on [a,b] × Rⁿ) of second order in the maximal grid size, which means that
  ||y_i-y(x_i)|| <= C × max{ x_i+1-x_i}²
  with an unknown constant C depending on y(x), F, R.
• We allow a variable grid size here, determined by a requirement on the local error like it is done in the step size control for initial value problems. We use here the technique of local extrapolation for this: the same problem is solved a second time, on a coherent grid with each subinterval didived in its middle. This results in a second solution y^{2}_i, i=1,...,2*N-1. Assuming second order convergence and an asymptotic development
  y_i = y(x_i)+ (h_i)²e₁(x_i) + O((h_i)⁴)
  (y^{2}_2i-1-y_i)/3 = err_i, i=1,...,N
  is an asymptotically correct estimate of the error in y_i.
  y^{2}_2i-1+err_i is the extrapolated solution which should be fourth order accurate.
  The user specifies two parameters tol and smally. Then we require that the weigthed error estimate
  yerror_i = ∑_{j=1 to n} |err_i,j|/max(|y^{2}_2i-1,j|,smally)
  satisfies
  yerror_i <= tol .
  If this is not the case, then the current subinterval i is divided into up to 8 subintervals for the next try. This is done in one sweep for all intervals. Thereby a new finer grid is obtained. The initial guess for the solution on this new finer grid is obtained interval by interval using quadratic interpolation of y^{2}, since we have there three grid values of this available in each interval.

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. These are all scalar second order problems of the type
  y'' = f(x,y,y'), g₁(y(a),y'(a)) = 0, g₂(y(b),y'(b)) = 0
  rewritten as first order systems of dimension 2 by substituting
  y₁(x) = y(x), y₂(x) = y'(x) .
  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.
• You select the number of grid points N by selecting N for the first (equidistant) grid. The grid size is then
    h = (b-a)/(N-1)
  
• Here n=2 and we automatically select the two components of y for plotting results. Since these components can have quite different sizes, each one is shown separately.
• 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 tol for the solution of this system and a parameter smally characterizing values of |y_i| which should not enter the relative error criterion (with the absolute precision requirement set to tol*smally) and printtable, which gives you the possibility to get a printed table of the solution on the last fine grid, or to set these values as desired. The default values are 30, 10^-4, 0.1, 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 final grid and the corresponding extrapolated value.
• The required computing time.
• If requested, a printed table of the computed solution on the last fine grid.
• (Up to) seven plots:
  1. the discretization errors on the last grid and for the extrapolated values for two selected components of y. In addition to the graph of this error you see (below the minimal function value) the generated grid (in green). For larger precision the grid however will be so fine that, due to the limited graphical presentation, there will appear a green bar.
  2. the development of the residual norm during the Newton process. The different curves correspond to the two runs and the different grids used in sequence red - green -blue etc
• Should the BVP-solver issue messages you get these as a protocol. For example if ||(h/2)*JF|| is larger than one, you will get a message that the grid should be refined. In addition, the refinement steps of the grid are reported. We limit the grid to 10000 points.

Questions: Predefined case?!

• Behaves the discretization error of the solver as expected ?
• In which cases you get a large pivotal growth and why?
• Which effect is to be seen if you take a too small N initially?
• What takes place if you require very small error (tol)?

to the input form: predefined case

Input: user defined case

• Here you have the possibility to define a BVP yourself. This requires
  1. Specification of n. We limit this by 10.
  2. Specification of a,b and of m_init with
       h = (b-a)/(m_init-1) for the first grid tried
     
  3. Whether to use some default parameters or not and whether to see a printed table of the solution.
  4. You write a piece of code for computing F and one more piece of code for computing R. 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 third 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 × m_init 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