Skip to content

Commit

Permalink
Rm scale in Laguerre and Hermite
Browse files Browse the repository at this point in the history
This is more elegantly handled by DiffMatOnDomain.
  • Loading branch information
amorison committed Feb 25, 2023
1 parent 90f89a1 commit a65875e
Show file tree
Hide file tree
Showing 4 changed files with 10 additions and 31 deletions.
33 changes: 6 additions & 27 deletions dmsuite/poly_diff.py
Original file line number Diff line number Diff line change
Expand Up @@ -244,34 +244,22 @@ class Hermite(DiffMatrices):
Attributes:
degree: Hermite polynomial degree, also the number of nodes.
scale: scaling factor.
"""

degree: int
scale: float # FIXME: this should be handled via composition

def __post_init__(self) -> None:
assert self.scale > 0.0

@property
def max_order(self) -> int:
"""Maximum order of derivative."""
return self.degree - 1

@cached_property
def norm_nodes(self) -> NDArray:
"""Hermite roots, unscaled."""
return herroots(self.degree)

@cached_property
def nodes(self) -> NDArray:
"""Scaled nodes."""
return self.norm_nodes / self.scale
return herroots(self.degree)

@cached_property
def _dmat(self) -> GeneralPoly:
# this is unscaled (scale == 1)
x = self.norm_nodes
x = self.nodes
alpha = np.exp(-(x**2) / 2) # compute Hermite weights.

# construct beta(l,j) = d^l/dx^l (alpha(x)/alpha'(x))|x=x_j recursively
Expand All @@ -285,7 +273,7 @@ def _dmat(self) -> GeneralPoly:
return GeneralPoly(at_nodes=x, weights=alpha, weight_derivs=beta[1:, :])

def at_order(self, order: int) -> NDArray:
return self.scale**order * self._dmat.at_order(order)
return self._dmat.at_order(order)


@dataclass(frozen=True)
Expand All @@ -296,33 +284,24 @@ class Laguerre(DiffMatrices):
Attributes:
degree: Laguerre polynomial degree. There are degree+1 nodes.
scale: scaling factor.
"""

degree: int
scale: float # FIXME: this should be handled via composition

@property
def max_order(self) -> int:
"""Maximum order of derivative."""
return self.degree

@cached_property
def norm_nodes(self) -> NDArray:
"""Laguerre roots, unscaled."""
def nodes(self) -> NDArray:
nodes = np.zeros(self.degree + 1)
nodes[1:] = lagroots(self.degree)
return nodes

@cached_property
def nodes(self) -> NDArray:
"""Scaled nodes."""
return self.norm_nodes / self.scale

@cached_property
def _dmat(self) -> GeneralPoly:
# this is unscaled (scale == 1)
x = self.norm_nodes
x = self.nodes
alpha = np.exp(-x / 2) # Laguerre weights

# construct beta(l,j) = d^l/dx^l (alpha(x)/alpha'(x))|x=x_j recursively
Expand All @@ -334,7 +313,7 @@ def _dmat(self) -> GeneralPoly:
return GeneralPoly(at_nodes=x, weights=alpha, weight_derivs=beta)

def at_order(self, order: int) -> NDArray:
return self.scale**order * self._dmat.at_order(order)
return self._dmat.at_order(order)


@dataclass(frozen=True)
Expand Down
4 changes: 2 additions & 2 deletions examples/laguerre.py
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
import matplotlib.pyplot as plt
import numpy as np

from dmsuite.poly_diff import Laguerre
from dmsuite.poly_diff import DiffMatOnDomain, Laguerre

lag = Laguerre(degree=31, scale=30.0)
lag = DiffMatOnDomain(xmin=0.0, xmax=3.0, dmat=Laguerre(degree=32))
x = lag.nodes
D1 = lag.at_order(1)
D2 = lag.at_order(2)
Expand Down
2 changes: 1 addition & 1 deletion tests/test_herdif.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
def test_herdif4() -> None:
"""Test of order 4 Hermite differentiation"""
expected = np.load("tests/data/herdif4.npy", allow_pickle=True)
herm = Hermite(degree=5, scale=1.0)
herm = Hermite(degree=5)
computed = np.zeros((3, herm.degree, herm.degree))
for order in range(1, 4):
computed[order - 1] = herm.at_order(order)
Expand Down
2 changes: 1 addition & 1 deletion tests/test_lagdif.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
def test_lagdif4() -> None:
"""Test of order 4 Laguerre differentiation"""
expected = np.load("tests/data/lagdif4.npy", allow_pickle=True)
lag = Laguerre(degree=4, scale=1.0)
lag = Laguerre(degree=4)
computed = np.zeros((3, lag.nodes.size, lag.nodes.size))
for order in range(1, 4):
computed[order - 1] = lag.at_order(order)
Expand Down

0 comments on commit a65875e

Please sign in to comment.