diff --git a/docs/volterra/volterra-example3.svg b/docs/volterra/volterra-example3.svg new file mode 100644 index 0000000..80247a6 --- /dev/null +++ b/docs/volterra/volterra-example3.svg @@ -0,0 +1,2054 @@ + + + + + + + + 2024-06-29T18:23:49.010682 + image/svg+xml + + + Matplotlib v3.9.0, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/example/volterra-example3.py b/example/volterra-example3.py new file mode 100644 index 0000000..ba0daf6 --- /dev/null +++ b/example/volterra-example3.py @@ -0,0 +1,46 @@ +import os + +import matplotlib.pyplot as plt +import numpy as np + +from inteq import SolveVolterra2 + +def k(s, t): + return 0.5 * (t-s)**2 * np.exp(s-t) +def free(t): + return -0.5 * t**2 * np.sin(t) * np.exp(-t/2) + +def true(t): + return 1/291762750 * np.exp( -3/2 * t ) * ( -3 * np.exp( t ) \ + * ( 16 * ( 2390928 + 365 * t * ( -9936 + 365 * t ) ) * np.cos( t ) \ + + 3 * ( 6965888 + 365 * t * ( 24544 + 25915 * t ) ) * np.sin( t ) \ + ) + 32 * ( 389017 * ( np.e )**( 3/2 * t ) + ( 3197375 * np.cos( \ + 1/2 * ( 3 )**( 1/2 ) * t ) + -318625 * ( 3 )**( 1/2 ) * np.sin( \ + 1/2 * ( 3 )**( 1/2 ) * t ) ) ) ) + +s, gmid = SolveVolterra2(k, free, a=0.0, b=20.0, num=100, method="midpoint") +s, gtrap = SolveVolterra2(k, free, a=0.0, b=20.0, num=100, method="trapezoid") +s, gsimp = SolveVolterra2(k, free, a=0.0, b=20.0, num=100, method="simpson") +s, ggreg = SolveVolterra2(k, free, a=0.0, b=20.0, num=100, + method="gregory", greg_order=7) + +figure, axs = plt.subplots(2, 1) +axs[0].plot(s, true(s)) + +axs[1].semilogy(s, np.abs(gmid-true(s))) +axs[1].semilogy(s, np.abs(gtrap-true(s))) +axs[1].semilogy(s, np.abs(gsimp-true(s))) +axs[1].semilogy(s, np.abs(ggreg-true(s))) + +axs[1].legend(["Midpoint", "Trapezoid", "Simpson", "Gregory(7)", "Exact"], loc="upper right") +axs[0].set_xlabel("t") +axs[0].set_ylabel("f(t)") +axs[1].set_xlabel("t") +axs[1].set_ylabel("Absolute error") + +figure.set_dpi(100) +figure.set_size_inches(8, 5) + +save_path = os.path.join(os.path.dirname(__file__), + "..", "docs", "volterra", "volterra-example3.svg") +figure.savefig(save_path, bbox_inches="tight") \ No newline at end of file diff --git a/inteq/quad.py b/inteq/quad.py new file mode 100644 index 0000000..cb652db --- /dev/null +++ b/inteq/quad.py @@ -0,0 +1,257 @@ +# SPDX-FileCopyrightText: 2024 Christopher Hillenbrand + +# SPDX-License-Identifier: MIT + +# MIT License + +# Copyright (c) 2024 Christopher Hillenbrand + +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: + +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. + +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + + +from fractions import Fraction +from functools import cache + +import numpy as np +import scipy.linalg as sla + + +@cache +def gregory_coef(k): + r"""Gregory coefficients. + + See https://en.wikipedia.org/wiki/Gregory_coefficients. + + Starting from :math:`G_0`, the first few coefficients are: + :math:`-\frac{1}{2}, \frac{1}{12}, -\frac{1}{24}, \frac{19}{720}, \ldots` + + .. math:: + G_k = \frac{-1}{(k+1)!} \left[ \frac{d^{k+1}}{dz^{k+1}}\frac{z}{\ln(1+z)}\right]_{z=0} + + Parameters + ---------- + k : int + (starts from zero) + + Returns + ------- + :class:`fraction.Fraction` + kth Gregory coefficient + """ + assert k >= 0 + if k == 0: + Fraction(-1, 2) + grn = Fraction(1, k + 2) + for r in range(k): + grn -= abs(gregory_coef(r)) * Fraction(1, k + 1 - r) + return grn * (-1) ** (k + 1) + + +@cache +def gregory_weights(k): + r"""Gregory quadrature weights. + + See equation 14 in + Fornberg, B., Reeger, J.A. An improved Gregory-like method for 1-D quadrature. + Numer. Math. 141, 1–19 (2019). https://doi-org.yale.idm.oclc.org/10.1007/s00211-018-0992-0 + + .. math:: + \mathbf{w} = + \begin{bmatrix} + 1 & 0 & 0 & 0 & \cdots \\ + -1 & 1 & 0 & 0 & \cdots \\ + 1 & -2 & 1 & 0 & \cdots \\ + -1 & 3 & -3 & 1 & \cdots \\ + \vdots & \vdots & \vdots & \vdots & \ddots \\ + \end{bmatrix} + \begin{bmatrix} + G_0 \\ + G_1 \\ + G_2 \\ + G_3 \\ + \vdots + \end{bmatrix} + + \begin{bmatrix} + 1 \\ + 1 \\ + 1 \\ + 1 \\ + \vdots + \end{bmatrix} + + Parameters + ---------- + k : int + Starts from 0; gives degree of polynomial approx. + + Returns + ------- + w : (k+1,) ndarray + Order k Gregory quadrature weights, exact. + """ + gr_coef = np.array([gregory_coef(r) for r in range(k + 1)]) + wts_exact = sla.invpascal(k + 1, kind="upper", exact=True) @ gr_coef + 1 + wts = np.asarray(wts_exact, dtype=np.float64) + wts.setflags(write=False) + return wts + + +@cache +def sigma_weights(k): + """Sigma quadrature weights (exact). + + See equation 98 in + + Schüler, M.; Golež, D.; Murakami, Y.; Bittner, N.; Herrmann, A.; + Strand, H. U. R.; Werner, P.; Eckstein, M. NESSi: + The Non-Equilibrium Systems Simulation Package. + Computer Physics Communications, 2020, 257, 107484. + https://doi.org/10.1016/j.cpc.2020.107484. + + See also + + P. H. M. Wolkenfelt, The Construction of Reducible Quadrature Rules for + Volterra Integral and Integro-differential Equations, IMA Journal of + Numerical Analysis, Volume 2, Issue 2, April 1982, Pages 131–152, + https://doi.org/10.1093/imanum/2.2.131. + + Parameters + ---------- + k : int + Starts from 0; gives degree of polynomial approx. + + Returns + ------- + sigma : (k+1, 2k+2) ndarray + Order k sigma quadrature weights, exact. + """ + + omega = gregory_weights(k) + sigma = np.zeros((k + 1, 2 * k + 2), dtype=object) + for i in range(1, k + 2): + sigma[i - 1, : k + 1] += omega + sigma[i - 1, i : k + 1] += np.flip(omega[i:]) - 1 + sigma[i - 1, k + 1 : k + 1 + i] = np.flip(omega[:i]) + sigma = np.asarray(sigma, dtype=np.float64) + sigma.setflags(write=False) + return sigma + + +@cache +def poly_integration_weights(k): + r"""Polynomial integration weights. + + .. math:: + I^{(k)} = \begin{bmatrix} + 0 & 0 & 0 & \cdots & 0 \\ + \frac{1}{1} & \frac{1}{2} & \frac{1}{3} & \cdots & \frac{1}{k+1} \\ + \frac{2^1}{1} & \frac{2^2}{2} & \frac{2^3}{3} & \cdots & \frac{2^{k+1}}{k+1} \\ + \vdots & \vdots & \vdots & \ddots & \vdots \\ + \frac{k^1}{1} & \frac{k^2}{2} & \frac{k^3}{3} & \cdots & \frac{k^{k+1}}{k+1} \\ + \end{bmatrix} (V^{(k)})^{-1} + + Parameters + ---------- + k : int + Polynomial degree. + + Returns + ------- + I : (k+1,k+1) ndarray + Polynomial integration weights. + """ + grid = np.arange(k + 1) + aux = np.vander(np.arange(k+1), increasing=True) / (grid+1) * grid[:, np.newaxis] + wts = aux @ np.linalg.inv(np.vander(np.arange(k+1), increasing=True)) + wts.setflags(write=False) + return wts + + +@cache +def gregory_weight_matrix(k, n): + """Weight matrix for Gregory quadrature. + + See equation 98 in + + Schüler, M.; Golež, D.; Murakami, Y.; Bittner, N.; Herrmann, A.; + Strand, H. U. R.; Werner, P.; Eckstein, M. NESSi: + The Non-Equilibrium Systems Simulation Package. + Computer Physics Communications, 2020, 257, 107484. + https://doi.org/10.1016/j.cpc.2020.107484. + + Parameters + ---------- + k : int + Polynomial degree. + n : int + Row (or upper limit of integration). + + Returns + ------- + w : (n+1,n+1) ndarray + Quadrature weights. + """ + nrow_max = max(n + 1, 2 * k + 2) + ncol_min = min(n + 1, 2 * k + 2) + assert n >= k + mat = np.zeros((nrow_max, n + 1), dtype=np.float64) + mat[:k + 1, :k + 1] = poly_integration_weights(k) + mat[k + 1 : 2 * k + 2, :ncol_min] = sigma_weights(k)[:, :ncol_min] + for r in range(2 * k + 2, n + 1): + mat[r, :r + 1] = 1 + mat[r, :k + 1] = gregory_weights(k) + mat[r, r - k : r + 1] = gregory_weights(k)[::-1] + mat = mat[:n + 1] + mat.setflags(write=False) + return mat + + +@cache +def simpson_quad(n): + """Repeated Simpson + Simpson's 3/8 rule + + Example 3.2.1 in + + H. Brunner & P. J. Van der Houwen, The Numerical Solution of + Volterra Equations, CWI Monographs, vol. 3, North-Holland, Amsterdam, 1986. + + Parameters + ---------- + n : int + Row (or upper limit of integration). + + Returns + ------- + w : (n+1,n+1) ndarray + Quadrature weights. + """ + arr = np.zeros((n + 1, n + 1), dtype=np.float64) + arr[1, 0] = 5/12 + arr[1, 1] = 8/12 + arr[1, 2] = -1/12 + array_3993 = np.array([3/8, 9/8, 9/8, 3/8], dtype=np.float64) + arr[2:, 0] = 1/3 + for row in range(2, n + 1, 2): + arr[row, 1:(2*row//2+1)] = np.tile([4/3, 2/3], row//2) + arr[row, row] = 1/3 + for row in range(3, n + 1, 2): + np.copyto(arr[row], arr[row - 3]) + arr[row, row-3:row+1] += array_3993 + arr.setflags(write=False) + return arr diff --git a/inteq/volterra.py b/inteq/volterra.py index eb6683d..f766737 100644 --- a/inteq/volterra.py +++ b/inteq/volterra.py @@ -6,6 +6,7 @@ import numpy import scipy.linalg +from .quad import simpson_quad, gregory_weight_matrix def solve( @@ -74,6 +75,7 @@ def solve2( b: float = 1.0, num: int = 1000, method: str = "midpoint", + greg_order: int = 3, ) -> numpy.ndarray: """ Approximate the solution, g(x), to the Volterra Integral Equation of the second kind: @@ -96,7 +98,9 @@ def solve2( num : int Number of estimation points between zero and `b`. method : string - Use either the 'midpoint' (default) or 'trapezoid' rule. + Quadrature method: 'midpoint' (default), 'trapezoid', 'simpson', or 'gregory'. + greg_order: int + Order of the Gregory quadrature rule. Used only if `method` is 'gregory'. Returns ------- @@ -109,22 +113,36 @@ def solve2( sgrid = numpy.linspace(a, b, num) h = (b - a) / (num - 1) # create a lower triangular matrix of kernel values - ktril = numpy.tril(k(sgrid, sgrid[:, numpy.newaxis])) + kmat = k(sgrid, sgrid[:, numpy.newaxis]) + fgrid = f(sgrid) if method == "midpoint": - pass + kmat = numpy.tril(kmat) + # first entry is exactly g(a) = f(a) + kmat[0, 0] = 0 + ggrid = scipy.linalg.solve_triangular( + numpy.eye(num) - h * kmat, fgrid, lower=True, check_finite=False + ) elif method == "trapezoid": # apply trapezoid rule by halving the left endpoints - ktril[:, 0] *= 1/2 + kmat = numpy.tril(kmat) + kmat[:, 0] *= 1/2 # apply trapezoid rule by halving the right endpoints (0,0 gets fixed later) - numpy.fill_diagonal(ktril, numpy.diag(ktril) / 2) + numpy.fill_diagonal(kmat, numpy.diag(kmat) / 2) + kmat[0, 0] = 0 + ggrid = scipy.linalg.solve_triangular( + numpy.eye(num) - h * kmat, fgrid, lower=True, check_finite=False + ) + elif method == 'simpson': + kmat *= simpson_quad(num - 1) + ggrid = scipy.linalg.solve(-h*kmat+numpy.eye(num), fgrid) + elif method == 'gregory': + kmat *= gregory_weight_matrix(greg_order, num - 1) + ggrid = scipy.linalg.solve(-h*kmat+numpy.eye(num), fgrid) else: - raise Exception("method must be one of 'midpoint', 'trapezoid'") - # first entry is exactly g(a) = f(a) - ktril[0, 0] = 0 + msg = f"Invalid method: {method}." + raise ValueError(msg) + # find the gvalues (/num) by solving the system of equations - ggrid = scipy.linalg.solve_triangular( - numpy.eye(num) - h * ktril, f(sgrid), lower=True, check_finite=False - ) # combine the s grid and the g grid return numpy.array([sgrid, ggrid])