Volterra Integral Equations¶
First Kind¶
This package provides the function inteq.SolveVolterra()
which approximates the solution to the Volterra Integral Equation of the first kind using the method in Betto and Thomas (2021).
-
inteq.
SolveVolterra
(k: Callable[[numpy.ndarray, numpy.ndarray], numpy.ndarray], f: Callable[[numpy.ndarray], numpy.ndarray] = <function <lambda>>, a: float = 0.0, b: float = 1.0, num: int = 1000, method: str = 'midpoint') → numpy.ndarray¶ Approximate the solution, g(x), to the Volterra Integral Equation of the first kind:
\[f(s) = \int_a^s K(s,y) g(y) dy\]using the method in Betto and Thomas (2021).
- Parameters
- kfunction
The kernel function that takes two arguments.
- ffunction
The left hand side (free) function with f(a) = 0.
- afloat
Lower bound of the integral, defaults to 0.
- bfloat
Upper bound of the estimate, defaults to 1.
- numint
Number of estimation points between zero and b.
- methodstring
Use either the ‘midpoint’ (default) or ‘trapezoid’ rule.
- Returns
- grid2-D array
Input values are in the first row and output values are in the second row.
Example¶
from inteq import SolveVolterra
import numpy as np
import matplotlib.pyplot as plt
# define kernel
def k(s, t):
return np.cos(s - t)
# define true solution
def trueg(s):
return (2 + s ** 2) / 2
# look decent with just 8 grid points but we'll use 25
s, g = SolveVolterra(k, num=25, method="trapezoid")
# plot functions
figure = plt.figure()
plt.plot(s, g)
plt.plot(s, trueg(s))
plt.legend(["Estimate", "True Value"])
plt.xlabel("s")
plt.ylabel("g(s)")
figure.set_dpi(100)
figure.set_size_inches(8, 5)
figure.savefig("..\\docs\\volterra\\volterra-example.svg", bbox_inches="tight")
Trapezoid vs Midpoint Rule¶
Volterra integral equations are typically solved using the midpoint rule. However, the trapezoid rule often converges faster. See below an example of the trapezoid rule performing well with just six grid points.
Thus, the trapezoid rule typically performs better. However, the trapezoid rule is less stable than the midpoint rule. An example where this this instability is an issue is provided below.
This can be remedied by smoothing the function. For example, with
inteq.helpers.smooth()
.
Second Kind¶
Implementation forthcoming.