Chasing

Introduction

The Chasing package provides a tool for numerical solution of linear multipoint boundary value problems. It has syntax of the built-in Mathematica function NDSolve and strongly enhances its functionality. In addition it provides a possibility to solve linear eigenfunction problem, when an initial guess range for eigenvalue is specified.

Description of the chasing method

Consider a system of boundary value problems in the form:

(Subscript[(∑_(i = 0)^N), k])(( (d^i) (Subscript[f, k]) (x))/(dx^i)) (Subscript[A, ik]) 
          (x) = 0        (k = 1, m); 
                          (1)

The linear boundary conditions are specified at some points inside the interval (xmin,xmax).
It is clear that this system can be rewritten in the form of system of first-order ODE with simple boundary conditions:

\(u' \((x)\)\  = A \((x)\)\ u \((x)\)\  + \ B \((x)\); 

    \(\(G\ \)\_i\) u \((x\_i)\)\  + \ b\_i\  = \ 0; \ 
    \((\(x\_1\  < \ x\_2\  < \ ...\) < x\_k)\).
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \((2)\)\)

Here u(x) denotes a vector function

u = {(Subscript[u, 1]), (Subscript[u, 2]), ..., (Subscript[u, N])}.

The main idea of the method is an assumption that for i-th boundary condition (specified at the point (Subscript[x, i]))  following linear relation is verified:

(Subscript[g, i]) (x) u (x)  +  (Subscript[h, i]) (x) = 0 , 
                                          (3)

where  (Subscript[g, i]) (x) is unknown vector function and (Subscript[h, i]) (x) is a free term. Differentiating equality (3) with respect to x yields:

(Subscript[g, i])' (x) u (x)  + (Subscript[g, i]) (x) u' (x)  +  
        (Subscript[h, i])' (x) = 
      0 .                                      (4)

Substitution of u'(x) from (2) leads to

\(g\_i' \((x)\) u \((x)\)\  + \(g\_i\) \((x)\) A \((x)\) u \((x)\)\  + 

        	\(g\_i\) \((x)\) B \((x)\) + \ h\_i' \((x)\) = 
      0\ .\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \((5)\)\)

Now we select one component of the vector u(x), say (Subscript[u, j]) , and express it through all other components of u, and functions g and h. ((Subscript[g, ij])(x) is set equal to 1, and its derivative to zero). The resulting expression we substitute into (4), and demand that the coefficients of the remaining components of u and the resulting free term are equal to zero. These conditions produce the auxillary problem to be solved in the forward step of integration. The number of equations in such auxillary problem is equal to dimension of the vector u - the dropped equation for (Subscript[g, ij])(x) is replaced by the equation for (Subscript[h, i])(x).
Each such problem is integrated from (Subscript[x, i]) to (Subscript[x, k]) = xmax. It results in a set of linear relations valid at the final point

(Subscript[g, i]) ((Subscript[x, k])) u ((Subscript[x, k]))  +  (Subscript[h, i]) ((Subscript[x, k])) = 
      0 .                               
        (6)

After adding to this set the original conditions at x=xmax, we arrive to a full set of conditions required to backward integration step. It should be noted here that the resulting set of conditions (linear algebraic system) should be resolved with respect to u ((Subscript[x, k])).  Sometimes this problem has no solution which corresponds to nonexistence of the solution for the original problem (1). Analysis of such cases is made using existing tools of Mathematica . 
	It is seen that this approach enables to deal with systems of linear boundary value problems (not only with a single problem). 
	A difficult problem may arise in some cases during forward integration procedures. Due to nonlinearity of the auxillary problem, it may appear, that the integration from (Subscript[x, i]) to xmax can not be made due to divergence of the solution inside the integration interval. There are two possibilities to resolve this difficulty. The first one is implemented in the current version of NDSolve - integration over some contour in complex plane, connecting (Subscript[x, i]) and xmax. 
	Another possibility is to use (less numerically accurate) linear version of chasing method. It is similar to the one described above, but there are no selection procedure for one component of vector u. It can be shown, that the relation (5) is verified if two following conditions are met (T denotes transposition of matrix):

\(g\_i' \((x)\)\  + \(A\^T\) \((x)\) \(g\_i\) \((x)\) = 0, \ 

    \(g\_i\) \((x)\) B \((x)\) + \ h\_i' \((x)\) = 0; 

    \(g\_i\) \((x\_i)\) = 
      \(\(G\_\(\(i; \)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\) \(h\_i\) \((x\_i)\) = 
        \(b\_i.\)\)\)

The approach depicted above can be also used for the numerical solution of eigenvalue problem. Consider a linear boundary value problem with additional parameter (eigenvalue) for which an interval of accessible values is specified. Then a kind of shooting method wrt to eigenvalue is organized which provides both the eigenvalue and eigenfunction for linear problem. This possibility also implemented in the current version of Chasing package.

Examples of usage

Loading the Package

At first add the directory containing it to the $Path variable (Note that MS Windows system needs a different directory separator)

AppendTo[$Path, "~/Here is the path to the required directory"];

and load the package.

Needs["Chasing`"]

Now you can look at short description of the only function Chasing.

?Chasing

Here is a set of options for Chasing.

Options[Chasing]

Simple Examples

Here is a typical example of the usage of the function Chasing. Note that its syntax is similar to that of NDSolve

\(sol3\  = \ 
      Chasing[{\(\(\(y'\)'\)'\)[x] + y[x] == 0, \(y'\)[0] == 0, 
          \(y'\)[2] == 1, \(\(y'\)'\)[0] == 0}, y[x], {x, 0, 2}]; 

    Plot[Evaluate[y[x] /. sol3], {x, 0, 2}, PlotRange -> All]; \)

You can check that the the solution is valid.

\(sol3a\  = \ 
      DSolve[{\(\(\(y'\)'\)'\)[x] + y[x] == 0, \(y'\)[0] == 0, 
          \(y'\)[2] == 1, \(\(y'\)'\)[0] == 0}, y[x], x]; 

    Plot[Chop[Evaluate[y[x] /. sol3a]], {x, 0, 2}, PlotRange -> All]; \)

Here are  more complicated problems are solved.

\(sol4a\  = \ 
      Chasing[{\(\(\(y'\)'\)'\)[x] + Sin[x]\ \(\(y'\)'\)[x] + x\ \(y'\)[x] == 
            0, y[0] == 1, \(y'\)[1] == 1, y[1] == 0}, y[x], {x, 0, 1}]; 

    Plot[Evaluate[y[x] /. sol4a], {x, 0, 1}, PlotRange -> All]; \)

\(sol4b\  = \ 
      Chasing[{\(\(\(y'\)'\)'\)[x] + Sin[x]\ \(\(y'\)'\)[x] + x\ \(y'\)[x] == 
            0, y[0] == 1, \(\(y'\)'\)[1] == 1, y[1] == 0}, y[x], {x, 0, 1}]; 
    
Plot[Evaluate[y[x] /. sol4b], {x, 0, 1}, PlotRange -> All]; \)

\(sol5\  = \ 
      Chasing[{\(\(y'\)'\)[x] + x\ \(y'\)[x] + y[x] == 0, y[0] == 1, 
          \(y'\)[1] == \(-0.7\)}, {y[x]}, {x, 0, 1}]; 

    Plot[Evaluate[y[x] /. sol5], {x, 0, 1}, PlotRange -> All]; \)

More Examples

The integration interval limits may not coinside with maximal and minimal points where boundary conditions are specified. Here two examples demonstrating the usage of Chasing in such cases.

\(sol6a\  = 
      Chasing[{\(\(y'\)'\)[x] + x\ \(y'\)[x] + y[x] == 0, y[0] == 1, 
          \(y'\)[1] == \(-0.7\)}, y[x], {x, 0, 0.5}]; 

    Plot[Evaluate[y[x] /. sol6a], {x, 0, 0.5}, PlotRange -> All]; \)

\(sol6b\  = 
      Chasing[{\(\(y'\)'\)[x] + x\ \(y'\)[x] + y[x] == 0, y[0] == 1, 
          \(y'\)[1] == \(-0.7\)}, y[x], {x, \(-1\), 0.5}]; 

    Plot[Evaluate[y[x] /. sol6b], {x, \(-1\), 0.5}, PlotRange -> All]; \)

Solution of several equations

NDSolve can solve a single boundary value problem. Chasing knows how to cope with a system of several such problems. Here is a simple example.

sol7  =  
      Chasing[{y''[x]  +  x z[x] + 5  y[x] == 0, 
          z'[x] == y[x] - 2  x, y'[0] == 0, y[1] == 1, 
          z[0.5] == 0.5}, {y[x], z[x]}, {x, 0, 1}]

Plot[Evaluate[{Re[y[x]], Re[z[x]]}  /.  sol7], {x, 0, 1}, 
      PlotRange -> All, PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}];

Solution of eigenvalue problem

Here we try to find an eigenfunction of the form y[x]=Sin[λ x] with λ=π. Note the low value for calculations tolerance specified by the option Tolerance. The form of output in this case slightly differs from that of NDSolve produces - an additional rule for eigenvalue is added.

sol8  =  
      Chasing[{y''[x] + λ^2 y[x] == 0, y[0] == 0, 
          y[1] == 0}, y[x], {x, 0, 1}, {λ, 2.8, 3.5}, 
        Tolerance  ->  4]

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

Example of nonstandard usage of Chasing - eigenfunction problem

If we call Chasing in standard way the results is similar to that of NDSolve produces, i.e., trivial solution.

\(sol1\  = \ 
      Chasing[{\(\(y'\)'\)[x] + 4  Pi^2\ y[x] == 0, y[0] == 0, y[1] == 0}, 
        y[x], {x, 0, 1}]; 

    Plot[Evaluate[y[x] /. sol1], {x, 0, 1}, PlotRange -> All]; \)

The reason of such behaviour that in this case we need to change the integration contour from the original to more complex in the complex plane. As a result the accuracy of the auxillary solution is very low, which produces the trivial solution. If you lower the tolerance of computation you may hope for nontrivial solution.

\(sol1\  = \ 
      Chasing[{\(\(y'\)'\)[x] + 4  Pi^2\ y[x] == 0, y[0] == 0, y[1] == 0}, 
        y[x], {x, 0, 1}, Tolerance -> 4]; 

    Plot[Evaluate[y[x] /. sol1], {x, 0, 1}, PlotRange -> All]; \)

You can also try to switch to another method of construction of auxillary system for slightly higher tolerance.

\(sol1\  = \ 
      Chasing[{\(\(y'\)'\)[x] + 4  Pi^2\ y[x] == 0, y[0] == 0, y[1] == 0}, 
        y[x], {x, 0, 1}, Tolerance -> 5, Method\  -> \ Linear]; 

    Plot[Evaluate[y[x] /. sol1], {x, 0, 1}, PlotRange -> All]; \)

Another example of nonstandard usage - problem with no solution

This example of a boundary value problem which has no solution at all. It can be checked with DSolve

DSolve[{y''[x] + Pi^2  y[x] == 0, y[0] == 1, y[1] == 0}, y[x], 
      x]

The solution produced by Chasing with standard calling is very bad (there should be no solution at all), but the reason is again the tolerance of calculations.

sol2 =  
      Chasing[{y''[x] + Pi^2  y[x] == 0, y[0] == 1, y[1] == 0}, 
        y[x], {x, 0, 1}]

If you lower the tolerance of computation you may hope for a proper result to be produced.

sol2 =  
      Chasing[{y''[x] + Pi^2  y[x] == 0, y[0] == 1, y[1] == 0}, 
        y[x], {x, 0, 1}, Tolerance -> 5]