We have tried to implement many of the examples found in Lund and Bower's book in SincLib, in order to measure the performance and correctness of SincLib. This page attempts to explain how to use various parts of the library. Once one feels comfortable programming specific instances of problems in SincLib, they might want to move on to implementing their own solvers.

We are always interested in seeing what kinds of problems the SincLib library is used for. Please contribute any problems you solve with SincLib back to the developers so we may include them in the standard set of examples.

Using Problem Files

In SincLib, the problem is written in a file separate from the solver implementation. For example, the DiffEq class solves the following problem:

We would like to solve this specific instance:

Create a file example_ivp.py with the following contents:

from numpy import *

def p( x ): return 0.0

def pp( x ): return 0.0

def q( x ): return x**-2.0

def f( x ): return (log(x) - 1.0)/x

def solution_func( x ): return x*log(x)

The first line imports the numpy math library and should be included in all problem files. The other functions defined here are:

• p: The function p as in the original ODE (set to 0 in this example).
  
• pp: The first derivative of p. While there are automatic differentiation utilities available in SincLib, it's usually a good idea to specify the derivative yourself. In this case, it is also 0.
  
• q: The function q.
  
• f: The right hand side of the ODE.
  
• solution_func: The solution to the ODE. While not necessary, this allows one to check the accuracy of the methods for known problems.

It is important to note that the variable passed to these functions is an array of values, not a single value. One should take care to make sure that any function used inside the function body should be one which is array-safe (most any math function should be).

The problem file describes all the functions and constants of the ODE. If a standard sinc constant is omitted (such as alpha, beta, or d), then the default is used. To solve this problem, simply type the following in a console after the SincLib environment has been set:

python DiffEq.py -p example_ivp.py

The -p option tells the class where to find the problem file (and should work for any class derived from SincBase. In fact, any class using the defaults system can override defaults using one of these problem files.). For a full list of command line parameters, use the -h flag. Here's the output of the online help:

$ python scripts/DiffEq.py -h
usage: DiffEq.py [options]

options:
  --homogeneous=HOMOGENEOUS
                        Boolean whether or not problem is homogenous
  --a=A                 Left boundary point.
  -u PP, --ppfunc=PP    
  --b=B                 Right boundary point.
  -d D, --d=D           The value for "d" in the sinc transform
  -f F, --function=F    
  -M M, --M=M           the number of nodes to use in the approximation
  -t TRANSFORM, --transform=TRANSFORM
                        The sinc transform to use.
  --b2=B2               Boundary condition at b
  -s SOLUTION_FUNC, --solution=SOLUTION_FUNC
                        The solution function for the problem.
  -r Q, --qfunc=Q       
  -q P, --pfunc=P       
  -b BETA, --beta=BETA  The value for beta in the sinc approximation
  --a2=A2               Boundary condition at a
  -p PROBLEM_FILE, --problemFile=PROBLEM_FILE
                        The file describing the second order Dirichlet-BC ODE
  -a ALPHA, --alpha=ALPHA
                        The value for alpha in the sinc approximation
  -h, --help            show this help message and exit

The output of the run using the problem file should be along these lines:

$ python scripts/DiffEq.py -p sinc.examples.example412
SincMethods.output:INFO:  Sinc error: 5.461636e-06

Then, the following picture should open (assuming X11 is running):

Basic Constants

The SincBase class, from which all our examples are derived, defines the following constants:

• M: the number of sinc nodes to use in the approximation.
• Default: 16
  
• beta: The value of beta from the sinc formulas.
• Default: 1
• alpha: The value of alpha from the sinc formulas.
• Default: 1
  
• transform: The conformal mapping to use in the problem.
• Default: eyelet
• d: The value for d in the sinc transform.
• Default: pi / 2

The constants alpha, beta, and d do not need to be touched if you don't quite understand what they are used for - basically, they affect the convergence properties of the sinc methods for specific cases.

The value of M controls how large the system is. There are (N+M+1) nodes in the system, where N = alpha / beta * M. Usually, the matrices involved in the computations have (N+M+1)^2D entries, where D is the dimension of the system (D=1 for ODEs, D=n for PDEs in dimension n, D=n for a system of n ODEs). In many cases, the convergence rate of the solution is exp(-M^k), where k=1/2 for quadrature, and make take other values for different problems. Refer to the literature. Note that direct solvers will have to keep all (N+M+1)^2D values of the matrix in memory; iterative solvers can compactly store the matrix, and may only need 2D(N+M+1) values in memory.

Since convergence and memory usage depend so heavily on M, it is usually better to not specify it in the problem file; instead, give it on the command line using the '-M' flag.

The conformal transform used depends heavily on the problem, and should usually be set in the problem file.

As each class explained below is a subclass from SincBase, all the above options are available, may be put in the problem file, and will not be repeated.

Quadrature in SincLib

The sinc.base.QuadratureLib module provides the Quadrature class, which solves definite integrals of the following form:
The problem file for the Quadrature class expects the following:

• integrand: The function which will be used for the integrand.
• answer: The exact answer to the integral, in case if you want to evaluate convergence rates.

The Quadrature class is the simplest example in SincLib. Notice that one does not define the limits of integration for this class; instead, these are defined implicitly by endpoints of the transform. For the identity transform (called none), the limits are -infinity to infinity. For the eyelet transform, the limits are 0 to 1.

If you just want a straightforward function which will perform the integration given the integrand and limits of integration (skipping over choosing the conformal mapping), try out the myQuad function found in the sinc.base.QuadratureLib module.

Example:

Save the following snippet of code to quad_example.py:

from numpy import *

M = 16                 # These two are not necessary, 
transform = 'eyelet'   # as they are the default setting

def integrand( x ): return 3*x**2

answer = 1.0

The problem outlined above is the integral of 3x² from 0 to 1 (implied by the use of the eyelet transform). The answer is, of course, 1.

Now, use the script Quadrature.py to solve the problem. The following snippet from the command line shows how the error sharply decreases from M=16 (the default) to M=64:

$ python scripts/Quadrature.py -p quad_example.py
Solution:  0.999993076187
Error:  6.92381261003e-06
$ python scripts/Quadrature.py -p quad_example.py -M 64 
Solution:  0.99999999997
Error:  2.97888380629e-11

Solving ODE IVPs

The sinc.ode module contains several submodules for solving ODEs. In the first section of this page, we showed an example using the DiffEq class, which solves a second-order ODE with Dirichlet boundary conditions. We now show several other available solvers.

Second-order ODE with Mixed Boundary Conditions

The module sinc.ode.DiffEqIc solves problems of the form:

It does this with the implemented class DiffEqIc. The following information is expected in the problem file:

• p, pp: The function p and its derivative, pp, in the original problem.
  
• q: The function q of the above problem.
  
• f: The right hand side of the problem.
  
• a, b: The boundary points of the ODE. These will not be necessary in later versions, as they are uniquely defined by the conformal transform.
  
• a0, a1, a2: The constants found in the boundary conditions for u(a).
  
• b0, b1, b2: The constants found in the boundary conditions for u(b).
  
• solution_func: The true solution of the problem, for convergence information.

Again, we defer to the Lund and Bower book for discussion of conditions placed upon the functions p, q, and f with respect to convergence rates. We also borrow an example from that book.

Example:

Save the following snippet of code to diffeqic_example.py:

from numpy import *

a = 0.0; b = 1.0
a0 = 1.0; a1 = 2.0; a2 = 1.0
b0 = 3.0; b1 = 1.0; b2 = 9.0

alpha = 1.0
beta = 1.0

def p( x ):
  return x**-1.0

def pp( x ):
  return -x**-2.0

def q( x ):
  return  x**-2.0

def f( x ):
  return sqrt(x) / 4.0 * ( -41.0*x**2.0 + 34.0*x - 1.0) - 2.0*x + x**-2.0

def solution_func( x ):
  return x**(5.0/2.0)*(1.0-x)**2 + x**3.0 + 1.0

Use the script DiffEqIc.py to run the solver:

$ python scripts/DiffEqIc.py -p diffeqic_example.py
SincMethods.output:INFO:  Sinc error: 3.900106e-04

It produces the following output:

First-order ODE with Mixed Boundary Conditions

The module sinc.ode.FirstOrderIcLib implements a first-order ODE solver for mixed boundary conditions. It solves equations of the form:
It expects the following information in the problem file:

• q: The function q from the first-order problem.
  
• f: The right-hand side of the first-order problem.
  
• solution_func: The true solution function; used to measure the accuracy of the sinc methods against known problems.

Example

The example for the first order IVP is:
The solution is:

Note that the value at infinity is 1, not 0. For this problem, we have the following file:

from numpy import *

transform='wedge'

def q(x): return -1.0

def f(x): return -1.0 + 2.0*exp(-x)

def solution_func(x): return 1.0 - exp(-x)

The script which handles the FirstOrderIc object is called FirstOrderIc.py. Here is an example session using the script:

$ python scripts/FirstOrderIc.py -p firstorderic_example.py 
SincMethods.output:INFO:  FirstOrderIc error at sinc nodes: 0.001505

M=N      h               AE              CR                      NOC 
16      0.785398163397  1.505016e-03    8.927596e-04    1.68580214827
32      0.55536036727   7.711328e-05    1.959845e-05    3.93466164564
64      0.392699081699  3.961882e-04    4.981374e-08    7953.39332176
128     0.277680183635  2.989736e-03    6.001552e-12    498160482.085

The output graph is the following:

Sturm-Liouville Problems

We now turn our attention to the Sturm-Liouville solver. We aim to solver problems of the form:

The output from this problem file is the following:

$ ./scripts/DiffEqOptv2.py -p sinc.examples.examples_easy
Solving using L-BFGS
Error of solution: 4.337839e-05

The graph is:

More Examples

The SincLib user walkthrough is continued on this page. It covers systems of ODEs and PDEs.