dedalus.libraries.dedalus_sphere.jacobi

Module Contents

dtype = 'float64'
coefficient_connection(N, ab, cd, init_ab=1, init_cd=1)

The connection matrix between any bases coefficients: Pab(z) = Pcd(z) @ Cab2cd –> Acd = Cab2cd @ Aab The output is always a dense matrix format. :Parameters: * N (int)

  • ab, cd (tuples of input and output Jacobi parameters.)

polynomials(n, a, b, z, init=None, Newton=False, normalised=True, dtype=dtype, internal='longdouble')

Jacobi polynomials, P(n,a,b,z), of type (a,b) up to degree n-1.

Jacobi.polynomials(n,a,b,z)[k] gives P(k,a,b,z).

Newton = True: cubic-converging update of P(n-1,a,b,z) = 0.

Parameters:

  • a,b > -1

  • z (float, np.ndarray.)

  • init (float, np.ndarray or None -> 1+0*z, or (1+0*z)/sqrt(mass) if normalised.)

  • normalised (classical or unit-integral normalisation.)

  • dtype (‘float64’,’longdouble’ output dtype.)

  • internal (internal dtype.)

quadrature(n, a, b, days=3, probability=False, dtype=dtype, internal='longdouble')

Jacobi ‘roots’ grid and weights; solutions to

P(n,a,b,z) = 0; len(z) = n.

sum(weights*polynomial) = integrate_(-1,+1)(polynomial), exactly up to degree(polynomial) = 2n - 1.

Parameters:

  • n (int > 0.)

  • a,b (float > -1.)

  • days (number of Newton updates.)

  • probability (sum(weights) = 1 or sum(weights) = mass(a,b))

  • dtype (‘float64’,’longdouble’ output dtype.)

  • internal (internal dtype (Newton uses by default).)

grid_guess(n, a, b, dtype='longdouble', quick=False)

Approximate solution to

P(n,a,b,z) = 0

measure(a, b, z, probability=True, log=False)

mu(a,b,z) = (1-z)**a (1+z)**b

if normalised: ((1-z)/2)**a ((1+z)/2)**b / (2*Beta(a+1,b+1))

Parameters:

a,b > -1

mass(a, b, log=False)

2**(a+b+1)*Beta(a+1,b+1) = integrate_(-1,+1)( (1-z)**a (1+z)**b )

Parameters:

  • a,b > -1

  • log (optional)

norm_ratio(dn, da, db, n, a, b, squared=False)

Ratio of classical Jacobi normalisation.

N(n,a,b) = integrate_(-1,+1)( (1-z)**a (1+z)**b P(n,a,b,z)**2 )

Gamma(n+a+1) * Gamma(n+b+1)

N(n,a,b) = 2**(a+b+1) * ———————————-

(2n+a+b+1) * Gamma(n+a+b+1) * n!

The function returns: sqrt(N(n+dn,a+da,b+db)/N(n,a,b))

This is used in rescaling the input and output of operators that increment n,a,b.

Parameters:

  • dn,da,db (int)

  • n (np.ndarray, int > 0)

  • a,b (float > -1)

  • squared (return N(n,a,b) or sqrt(N(n,a,b)) (defalut))

operator(name, normalised=True, dtype=dtype)

Interface to base JacobiOperator class.

Parameters:

  • name (A, B, C, D, Id, Pi, N, Z (Jacobi matrix))

  • normalised (True –> unit-integral, False –> classical.)

  • dtype (output dtype)

class JacobiOperator(name, normalised=True, dtype=dtype)

The base class for primary operators acting on finite row vectors of Jacobi polynomials.

<n,a,b,z| = [P(0,a,b,z),P(1,a,b,z),…,P(n-1,a,b,z)]

P(k,a,b,z) = <n,a,b,z|k> if k < n else 0.

Each oparator takes the form:

L(a,b,z,d/dz) <n,a,b,z| = <n+dn,a+da,b+db,z| R(n,a,b)

The Left action is a z-differential operator. The Right action is a matrix with n+dn rows and n columns.

The Right action is encoded with an “infinite_csr” sparse matrix object. The parameter increments are encoded with a JacobiCodomain object.

L(a,b,z,d/dz) ………………………. dn, da, db

A(0) = a ………………………. 0, 0, 0

B(+1) = 1 ………………………. 0, 0, +1 B(-1) = 1+z ………………………. +1, 0, -1 B(0) = b ………………………. 0, 0, 0

C(+1) = b + (1+z)d/dz ………………… 0, +1, -1 C(-1) = a - (1-z)d/dz ………………… 0, -1, +1

D(+1) = d/dz ……………………….. -1, +1, +1 D(-1) = a(1+z) - b(1-z) - (1-z)(1+z)d/dz .. +1, -1, -1

Each -1 operator is the adjoint of the coresponding +1 operator.

In addition there are a few exceptional operators:

Identity: <n,a,b,z| -> <n,a,b,z|

Parity: <n,a,b,z| -> <n,a,b,-z| = <n,b,a,z| Pi

The codomain is not additive in this case.

Number: <n,a,b,z| -> [0*P(0,a,b,z),1*P(1,a,b,z),…,(n-1)*P(n-1,a,b,z)]

This operator doesn’t have a local differential Left action.

Variables:

  • name (str) – A, B, C, D

  • normalised (bool) – True gives operators on unit-integral polynomials, False on classical normalisation.

  • dtype ('float64','longdouble') – output dtype.

__call__(p): p=-1,0,1

returns Operator object depending on p. Operator.function is an infinite_csr matrix constructor for n,a,b. Operator.codomain is a JacobiCodomain object.

staticmethods()
-------------
identity: Operator object for identity matrix
parity:   Operator object for reflection transformation.
number:   Operator object for polynomial degree.
property normalised
dtype = 'longdouble'
static identity(dtype=dtype)
static parity(dtype=dtype)
static number(dtype=dtype)
class JacobiCodomain(dn=0, da=0, db=0, pi=0, Output=None)

Base class for Jacobi codomain.

codomain = JacobiCodomain(dn,da,db,pi)

n’, a’, b’ = codomain(n,a,b)

if pi == 0:

n’, a’, b’ = n+dn, a+da, b+db

if pi == 1:

n’, a’, b’ = n+dn, b+da, a+db

pi_0 + pi_1 = pi_0 XOR pi_1

Variables:

__arrow (stores dn,da,db,pi.)

self[0:3]: returns dn,da,db,pi respectively.
str(self): displays codomain mapping.
self + other: combines codomains.
self(n,a,b): evaluates current codomain.
-self: inverse codomain.
n*self: iterated codomain addition.
self == other: determines equivalent codomains (a,b,pi).
self | other: determines codomain compatiblity and returns larger-n space.