Solving Systems of ODEs

The ODE system solver is based on ODE-IVP-PACK by Stenger et al. The algorithm is presented in this paper:

• F. Stenger, B. Keys, M. O'Reilly, K.Parker, and S.-A Gustafson. OPE-IVP-PACK via Sinc indefinite integration and Newton's method NUMERICAL ALGORITHMS, 20 (1999), PP. 241-268

The sinc.base.ode_system module provides many ODE solver submodules. For the system of ODEs, the best module to use is DiffEqSystemv2Lib; this provides the DiffEqSystem2 class. It solves the following system:
The inputs needed are:

• equations: The number of equations, n, of the system.
  
• init: The initial conditions for the equations; should be an n-tuple of floating point values.
  
• yinit: The initial guess for the solution function.
  
• f: The right hand side of the equation.
  
• j: The Jacobian of f.
  
• solution_func: The true solution.
  

Example

Below is a fully-commented problem-file for a system of ODEs.

# This problem is the simple second order IVP:

# y'' + y = 0, y(0) = 0

# The solution is:

# y(x) = sin(x) 

# Converted to a system of ODEs, we have:

# y1' = y2
# y2' = -y1
# y1( 0 ) = 0; y2( 0 ) = 1

# Where the solution is:

# y1(t) = sin(t)
# y2(t) = cos(t)

# The standard line to import the numerical library

from numpy import *

equations = 2

init = (0.0,1.0)

def yinit( t ):
  # Set the initial solution guess to be two constant functions,
  # y1 := 1 and y2 := 0
  y1 = 0.0
  y2 = 1.0
  return array(( y1, y2 ))

def f( t, y ):
  y1 = y[1]
  y2 = - y[0]
  return array(( y1, y2 ))

def j( t, y ):
  Bij = zeros((2,2,t.shape[0]), float64 )
  Bij[0,0] = 0.0
  Bij[0,1] = 1.0
  Bij[1,0] = -1.0
  Bij[1,1] = 0.0
  return Bij

def solution_func( t ):
  y1 = sin( t )
  y2 = cos( t )
  return array((y1, y2))

The command line needed to execute this script (Also found in the SincLib distribution as the module sinc.examples.sine_sys):

$ ./scripts/DiffEqSystem.py -p sinc.examples.sine_sys -M 32
Total time: 0.470304965973
Compute time: 0.456157922745
Prep time: 0.0141470432281
Relative Sinc error in our solution: 5.724213e-08 

The resulting graph is:

Solving PDEs

The sinc.pde module contains submodules which solve various specific PDEs. This section shows examples of solving the Heat equation, Poisson's equation, and Schrodinger's equation.

Heat Equation

The heat equation that we are interested in is:


We implement this in the sinc.pde.HeatLib with the Heat class. The following inputs are required for the input file:

• dimensions: The dimensions of the problem
  
• f: The right hand side for the heat equation.
  
• solution_func: The true solution to the problem; used for measuring the convergence rate of the solver.
  

Example

The following problem file is included as the module sinc.examples.example62:

from scipy import *

dimensions = 2

alpha = (.5,.5)
beta = (.5, .5)
M = 4

def f( x, t ):
  expt = exp(-t)
  xterm = (x*(1.0-x)*(1.0-t) + 2.0*t)
  return expt * xterm

def solution_func( x, t ):
  return x*(1.0-x)*t*exp(-t)

The heat equation example is run by the Heat.py script. The output is as follows:

$ python ./scripts/Heat.py -p sinc.examples.example62 -M 8
SincMethods.Heat:INFO:  Starting compute_iterative for the Heat problem.
SincMethods.Heat:DEBUG:  Starting compute
SincMethods.Heat:DEBUG:  About to make nnodes.
SincMethods.Heat:DEBUG:  Made nodes.
SincMethods.Heat:DEBUG:  Applied f.
SincMethods.Heat:DEBUG:  Made kron; shape: (289, 289)
SincMethods.Heat:DEBUG:  Finished with left-hand-side.
SincMethods.Heat:DEBUG:  Final problem size: (289, 289)
SincMethods.Heat:INFO:  Starting solve (time for setup: 0.069798)
SincMethods.Heat:DEBUG:  Return value from GMRES method is 0
SincMethods.Heat:DEBUG:  Residual from GMRES is: 2.237370e-07
SincMethods.Heat:INFO:  Finished solution in 0.443027 seconds
SincMethods.Heat:DEBUG:  Max difference: Uniform, 0.007486; Sinc, 0.001821
Hit enter to continue.

After the last line, if gnuplot is installed on the computer, the following 3d image appears:

If gnuplot 4.0 or higher is installed, you will be able to grab and turn the image with your mouse.

Poisson's Equation

The equation we are interested in here is:

We implement this in the sinc.pde.PoissonLib with the Poisson class. The inputs required are:

• dimensions: Number of dimensions for the problem.
  
• f: The left hand side of the equation.
  
• solution_func: The true solution.
  

Example

Below is an example problem file for the poisson solver:

from scipy import *
import logging

dimensions = 2

alpha = (.5,.5)
beta = (.5, .5)
M = 2

def f( x, y ):
  a = log(y)
  a /= x
  a *= -y
  b = - x*log(x)/y
  return a+b

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

The above problem is distributed as sinc.examples.example56, and it is example 5.6 from Lund and Bower's book.
The poisson solver is run with the script Poisson.py in the scripts/ directory. Here's an example command line session:

$ ./scripts/Poisson.py -p sinc.examples.example56 -M 6   
root:INFO PoissonLib:61:  Finished common part in 0.002722 seconds
SincMethods.Poisson:INFO:  Starting compute_iterative.
SincMethods.Poisson:INFO:  Starting solve (time for setup: 0.004402)
root:INFO PoissonLib:105:  Finished solution in 0.060672 seconds
Hit enter to continue.

The resulting graph is:

Schrodinger's equation

We are interested in solving the time-independent Schrodinger's equation:


We implement this in the module sinc.pde.SchrodingerLib with the class Schrodinger. The following inputs are expected:

• dimensions: The dimension of the problem.
  
• V: The potential function, V, in the equation.
  
• base_eigenvalue: The smallest eigenvalue of the system.
• base_eigenfunc: The eigenfunction associated to the base_eigenvalue.
  

Example

from scipy import *

dimensions = 2

alpha = (.5,.5)
beta = (.5, .5)
M = 6

# The constraints on the potential are
# put into place by the fact we are working
# on the interval 0 <= x <= 1
def v( x, t ):
  return 0*x

base_eigenvalue = pi**2*2

def base_eigenfunc( x, y ):
  return 2*sin(pi*x)*sin(pi*y)

The above example is available as sinc.examples.schrodinger_box. The Schrodinger equation solver is Schrodinger.py; here is a sample run:

$ python ./scripts/Schrodinger.py -p sinc.examples.schrodinger_box
root:INFO SchrodingerLib:66:  Finished common part in 0.004692 seconds
SincMethods.Schrodinger:INFO:  Starting compute_direct.
SincMethods.Schrodinger:INFO:  Final problem size: (169, 169)
SincMethods.Schrodinger:INFO:  Starting solve (time for setup: 0.052832)
SincMethods.Schrodinger:INFO:  Finished solution in 0.933610 seconds
148 (19.3521981747+0j)
SincMethods.Schrodinger:INFO:  Relative inf error: 0.059898
SincMethods.Schrodinger:INFO:  Difference between actual and computed base energies: -3.870106e-01

The output graph is:

Solving PDEs with the Sinc Method of Lines

A method of lines-based solver is found in the module sinc.pde.SincMolLib, as the class SincMol. The method of lines can be done in three different ways:

• Fully-Sinc approach; the temporal and spatial dimensions are modeled using sinc methods.
• Partial-Sinc approach 1; the temporal dimension is solved using an industry-standard solver, but the spatial dimension is modeled with sinc methods.
• Partial-Sinc approach 2; the spatial coordinates are modeled with forward-differences, but the temporal dimension is solved with sinc methods.

Currently, SincMol only implements the first approach (however, it should be noted that implementations of the second two approaches are available, but not yet satisfactory). When sinc methods are used for the temporal dimension, we model the PDE as a system of linear ODEs using the methods outlined above. Currently, only first-order PDEs are supported (i.e., the transport equations). The generic problem the SincMol is:
The example problem which we solve is:
The resulting problem file, which is distributed as sinc.examples.first_example, is:

from numpy import *

def f( x ): return exp(-x**2)

def p( t, x ): return 1.0

def q( t, x ): return 2.0

def solution_func( t, x ):
  return exp(-(x-t)**2-2*t)

The driver script in the scripts/ directory is called SincMol.py. A sample command-line session is show below:

The resulting graph is:
Each line on the graph represents a certain timeslice; the actual and computed values are shown (but the lines are nearly-overlapping). You can see that the solution graph decreases exponentially in time, as expected from the solution graph.