solve the Fredholm integral equation of the first kind

from inteq import SolveFredholm
import numpy as np
import matplotlib.pyplot as plt

miu = 0.2

# define kernel
def k(s, t):
    return 1/(1+100*(t-s)**2)


# define free function
def f(s):
    return np.exp(-(s-0.5)**2/0.06)


# define true solution
def trueg(s):
    return k(0, s)


# apply the solver
s, g = SolveFredholm(k, f, a=-1, b=1, num=1000, smin=-1, smax=1)

# 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\\fredholm\\fredholm-example.svg", bbox_inches="tight")


plt.show()

# -*- coding: utf-8 -*-

import typing
import numpy
from .helpers import makeH, simpson

#%%


def solve(
    k: typing.Callable[[numpy.ndarray, numpy.ndarray], numpy.ndarray],
    f: typing.Callable[[numpy.ndarray], numpy.ndarray] = lambda x: x,
    a: float = -1.0,
    b: float = 1.0,
    gamma: float = 1e-3,
    num: int = 40,
    **kwargs,
) -> numpy.ndarray:
    """
    Approximate the solution, g(x), to the Fredholm Integral Equation of the first kind:
    
    $$
    f(s) = \\int_0^b K(s,y) g(y) dy
    $$

    using the method described in Twomey (1963). It will return a smooth curve that is an approximate solution. However, it may not be a good approximate to the true solution.

    Parameters
    ----------
    k : function
        The kernel function that takes two arguments.
    f : function 
        The left hand side (free) function that takes one argument.
    a : float
        Lower bound of the of the Fredholm definite integral, defaults to -1.
    b : float
        Upper bound of the of the Fredholm definite integral, defaults to 1.
    num : int
        Number of estimation points between zero and `b`.
    smin : float
        Optional. Lower bound of enforcement values for s.
    smax : float
        Optional. Upper bound of enforcement values for s.
    snum : int
        Optional. Number of enforcement points for s.

    Returns
    -------
    grid : 2-D array
        Input values are in the first row and output values are in the second row.
    """
    if not isinstance(num, int):
        num = int(num)
    # need num to be odd to apply Simpson's rule
    if (num % 2) == 0:
        num += 1
    # set defaults for s params
    if "smin" in kwargs.keys():
        smin = kwargs["smin"]
    else:
        smin = a
    if "smax" in kwargs.keys():
        smax = kwargs["smax"]
    else:
        smax = b
    if "snum" in kwargs.keys():
        snum = kwargs["snum"]
    else:
        snum = 2 * num
    # get grid for s
    sgrid = numpy.linspace(smin, smax, snum)
    # get grid for y
    ygrid = numpy.linspace(a, b, num)
    # get quadrature weights
    weights = simpson(num)
    # create the quadrature matrix
    ksqur = weights * k(sgrid[:, numpy.newaxis], ygrid)
    # Make H matrix as in (Twomey 1963)
    if gamma != 0:
        Hmat = makeH(num)
        AAgH = (ksqur.T @ ksqur) + (gamma * Hmat)
    else:
        AAgH = ksqur.T @ ksqur
    # find the gvalues (/num) by solving the system of equations
    ggrid = numpy.linalg.solve(AAgH, ksqur.T @ f(sgrid))
    # combine the s grid and the g grid
    return numpy.array([ygrid, (ggrid * num / (b - a))])


# %%

import numpy

#%%
def makeH(dim: int) -> numpy.ndarray:
    """
    Make H matrix for estimating Fredholm equations as in (Twomey 1963).

    Parameters
    ----------
    dim : int
        The dimension of the H matrix.

    Returns
    -------
    Hmat : 2-D array
        The H matrix 
    """
    # create (symmetric) diagonal vectors
    d2 = numpy.ones(dim - 2)
    d1 = numpy.concatenate(([-2], numpy.repeat(-4, dim - 3), [-2]))
    d0 = numpy.concatenate(([1, 5], numpy.repeat(6, dim - 4), [5, 1]))
    # make matrix with diagonals
    Hmat = numpy.diag(d1, 1) + numpy.diag(d2, 2)
    Hmat = Hmat + Hmat.T
    numpy.fill_diagonal(Hmat, d0)
    return Hmat


#%%
def simpson(dim: int) -> numpy.ndarray:
    """
    Make H matrix for estimating Fredholm equations as in (Twomey 1963).

    Parameters
    ----------
    dim : int
        The dimension of the H matrix.

    Returns
    -------
    weights : 1-D array
        The quadrature weights according to Simpson's rule.
    """
    if dim > 2 and (dim % 2) == 1:
        tiles = numpy.tile([4, 2], dim // 2)
        return numpy.concatenate(([1], tiles[0 : dim - 2], [1])) / 3
    else:
        raise ValueError("Simpson's rule requires an odd number of endpoints")


#%%
def smooth(v: numpy.ndarray) -> numpy.ndarray:
    """
    Smooth a vector that is oscillating.

    Parameters
    ----------
    v : 1-D array
        The oscillating vector you want to smooth.

    Returns
    -------
    sv : 1-D array
        The smoothed vector.
    """
    dim = len(v)
    smat = numpy.diag(
        numpy.concatenate(([1], numpy.repeat(0.5, dim - 1)))
    ) + numpy.diag(numpy.repeat(0.5, dim - 1), -1)
    return smat @ v