Theory

The ODE we are interested in solving is:

The resulting Sinc-Galerkin System is:

Each of the terms in the above matrix system can be expressed as combinations of the matrices I0, I1, or I2 provided by the sinclib module or the function f1, f2, f3 provided by the defaults system.  The convergence theorem provided from Lund and Bower's book guarantees the following convergence rate:

Implementation

There are a few important parts to computing the above system:

• Pre-computing all the diagonal functions, such as f1, f2, f3, and the Sinc matrices I0, I1, I2.
  
• Computing the matrix A.
• Assemble the discrete system.
• Solve the matrix equations.

	def update(self):
		""" Update the basic properties of the system. """
		super( DiffEq, self ).update()
		M = self.M; N = self.N
		self.I0 = cSincLib.genSincMat(M+N+1,0)
		self.I1 = cSincLib.genSincMat(M+N+1,1)
		self.I2 = cSincLib.genSincMat(M+N+1,2)
		self.df1, self.df2, self.df3 = None, None, None
		self.updateCoeff()
		self.nodes = cSincLib.genSincNodes( self.psi, self.M, self.N, self.h )
		self.df1 = diag(map(self.f1,self.nodes))
		self.df2 = diag(map(self.f2,self.nodes))
		self.df3 = diag(map(self.f3,self.nodes,[1.0]*len(self.nodes)))
		self.sA = (-1.0/self.h**2)*self.I2 \
			+ (1.0/self.h)*dot(self.I1, self.df1) + self.df3

The next function, compute_base, pre-computes the vectors for functions q, p, and the right hand side of the discrete system (saved to the variable b).

	def compute_base( self ):
		""" The common base of all the newer computation methods. """
		self.p2nodes = self.pp( self.nodes )*self.f2( self.nodes )  \
			- self.p( self.nodes )*self.f1( self.nodes )
		self.qnodes = self.q( self.nodes )
		self.pnodes = self.p( self.nodes )
		self.b = self.df2*self.df2*self.df

After compute_base is run, the function compute_new is called for direct solves (compute_iterative_new is called if iterative methods are requested):

	def compute_new( self ):
		""" A better method for computing the solution using C-speedups.
		It still uses direct methods, however - not good for large matrices.
		"""
		self.compute_base()
		self.tmp2 = -self.df2 * self.p2nodes + numpy.power(self.df2,2.0) \
			* self.qnodes
		self.tmp1 = self.pnodes * self.df2
                self.part1 = -1/self.h*cSincLib.mdMultiply(self.I1,self.tmp1)
                self.parts = cSincLib.mdAdd( self.part1, self.tmp2)
		self.A = self.sA +  self.parts
		self.soln_vector = solve(self.A, self.b)
		self.solution = self.make_solution( self.soln_vector )
		return self.solution

compute_new assembles the discrete Sinc-Galerkin system.  The three lines calls solve, which solves the system Au = b, then creates the solution function with the method self.make_solution.  This takes the solution vector, self.soln_vector, and creates the following solution function:

Here, the elements of the solution vector are the constants u_k.  The solution function is then returned.