boundary_value_problem_1D_engine.py
No OneTemporary
Actions

Subscribers

None

File Metadata

	import numpy as np
	import scipy.sparse as scsp
	from rrompy.hfengines.scipy_engines import EigenproblemEngine

	class BVP1DEngine(EigenproblemEngine):
	"""
	Second order finite differences (uniform grid) for
	y''(x)+ky(x)=f(x) (0<x<1)
	y(0)=0
	y'(0)+ay'(1)=0 (0<a<1)
	((( //math.nist.gov/MatrixMarket/data/NEP/mvmode/mvmode.html )))
	"""
	def __init__(self, a:float, n:int, shift : np.complex = 0.,
	seed : int = 42):
	h = 1. / (n + 1)
	A0 = scsp.diags([-2. * np.ones(n - 1), np.ones(n - 1), np.ones(n - 2)],
	[0, 1, -1], shape = (n - 1, n),
	format = "csr", dtype = np.complex)
	a0l = np.zeros((1, n))
	a0l[0, : 2] = [4., -1.]
	a0l[0, -3 :] = [a, -4. * a, 3. * a]
	A0 = scsp.vstack([A0, scsp.csr_matrix(a0l, dtype = np.complex)])
	A1 = scsp.diags([np.append(h ** 2. * np.ones(n - 1), 0.)], [0],
	shape = (n, n), format = "csr", dtype = np.complex)
	super().__init__([A0 + shift * A1, A1], seed)

	def BVP1DPoles(a:float, km:float, kM:float, shift : np.complex = 0.):
	N = np.ceil(.5 * (kM ** .5 / np.pi - 1.))
	sigma = 1. / a + (1. / a 2. - 1) .5
	k = np.arange(- N, N)
	pls = ((((2 * k + 1) * np.pi) 2. - np.log(sigma) 2.)
	+ 2.j * (2 * k + 1) * np.pi * np.log(sigma)) - shift
	return pls[np.logical_and(np.real(pls) >= km, np.real(pls) <= kM)]