dymad.numerics.weak

Functions

compute_newton_cotes_weights(num_points, dx, ...)	Compute Newton-Cotes integration weights for a set of num_points samples.
generate_weak_weights(dt, n_integration_points)	Generate weights for a weak formulation approach, returning arrays (C, D) and the number of "windows" K along the time dimension.
jacobi_polynomial(order, coords)	Evaluate the Jacobi polynomials of order order (from 0 to order-1) with parameters alpha = 1, beta = 1 at the given coordinates.
jacobi_polynomial_derivative(order, coords)	Evaluate the first derivative of the Jacobi polynomials with alpha=1, beta=1 for orders 1 to order-1, at the given coordinates.

dymad.numerics.weak.compute_newton_cotes_weights(num_points, dx, order)

Compute Newton-Cotes integration weights for a set of num_points samples.

This function assumes that (num_points - 1) is divisible by order. For example, for Simpson’s 1/3 rule (order=2), we need (num_points - 1) to be even.

Parameters:

• num_points (int) – Total number of points (N).

• dx (float) – The spacing between consecutive sample points.

• order (int) – Newton-Cotes rule to use (1, 2, 3, or 4).

Returns:

A 1D array of length num_points containing the integration weights.

Return type:

np.ndarray

dymad.numerics.weak.generate_weak_weights(dt, n_integration_points, poly_order=4, int_rule_order=4)

Generate weights for a weak formulation approach, returning arrays (C, D) and the number of “windows” K along the time dimension.

Steps:

1. Compute a length scale:

\[L = \frac{t_{N-1} - t_0}{2}\]

2. Generate Jacobi polynomial basis (P0) and its derivative (P1), each with alpha=1, beta=1, on a grid of size n_integration_points in [-1,1].

3. Construct weighting arrays for the integrals:

\[w_0 = 1 - h^2,\quad w_1 = -\frac{2h}{L}\]

4. Compute Newton-Cotes integration weights on that grid.

5. Combine everything to form:

\[C = -(P_1 w_0 + P_0 w_1) w,\quad D = P_0 w_0 w\]

Parameters:

• dt (float) – Time step size.

• n_integration_points (int) – Number of integration points (N) in [-1, 1].

• poly_order (int) – Order of the Jacobi polynomial basis.

• int_rule_order (int) – Order of Newton-Cotes integration rule (1..4).

Returns:

Weights C and D:

• C: shape (poly_order, n_integration_points),

• D: shape (poly_order, n_integration_points).

Return type:

Tuple[np.ndarray, np.ndarray]

dymad.numerics.weak.jacobi_polynomial(order, coords)

Evaluate the Jacobi polynomials of order order (from 0 to order-1) with parameters alpha = 1, beta = 1 at the given coordinates.

Parameters:

• order (int) – The highest polynomial order to evaluate (exclusive). We compute polynomials 0, 1, …, order-1.

• coords (np.ndarray) – 1D array of points at which to evaluate.

Returns:

A 2D array of shape (order, len(coords)) where each row is the evaluated polynomial of a specific order.

Return type:

np.ndarray

dymad.numerics.weak.jacobi_polynomial_derivative(order, coords)

Evaluate the first derivative of the Jacobi polynomials with alpha=1, beta=1 for orders 1 to order-1, at the given coordinates.

Note

The returned array has shape (order, len(coords)) for consistency, but the 0-th row remains all zeros (since derivative is not computed for n=0).

Parameters:

• order (int) – The highest polynomial order to evaluate (exclusive).

• coords (np.ndarray) – 1D array of points at which to evaluate.

Returns:

A 2D array of shape (order, len(coords)) where row i holds the derivative of Jacobi polynomial i at those coordinates. Row 0 is zeros by definition here.

Return type:

np.ndarray