How to program the convergence of a sequence of systems of integral equations using Scipy

Question

I'm trying to solve the problem

where $u_n$ and $v_n$ are sequences that converge to the solution $u$ and $v$ and $\lambda$, $\sigma$, $f$ and $g$ and K_lambda are all given.

I thought of using the scipy module integrate with a while loop, but obviously I don't know how to implement the condition of convergence. I also don't know how to implement u_n(t) as integrand. Here's what I tried.

import numpy as np
import scipy.integrate as integrate
import scipy.special as special
import matplotlib
import matplotlib.pyplot as plt
import math

N = 30


u0 = np.ones(N)*0.2
v0 = np.ones(N)*0.15

bta = 2.1
eps = 0.1
sgm = 1.5

def K(t, s, lbda, T):
    L1= np.exp(lbda*(T-s))*np.exp(lbda*t)/(1-np.exp(lbda*T))
    if t>=s:
      L2 = np.exp(lbda*(t-s))
    else :
      L2=0.0
    return L1+L2


def f(s,u,v):
  return np.sin(s)**2-bta*(u/(1+v))
def g(s,u,v):
  return -np.cos(s)**2+sgm*(u/(1+v))+eps/(1+v)

def sigma(s):
    return 1+np.sin(s)**2

def integrand_u(s, t, lbda, T,u,v):
    return K(t, s, lbda, T) * (sigma(s) + lbda*u + f(s,u,v))

def integrand_v(s, t, lbda, T,u,v):
    return K(t, s, lbda, T) * (sigma(s) + lbda*v + g(s,u,v))

def integral_u(t, lbda, T,u,v):
    result = integrate.quad(integrand_u, 0, T, args=(t, lbda, T, u, v))[0]
    return result

def integral_v(t, lbda, T,u,v):
    result = integrate.quad(integrand_v, 0, T, args=(t, lbda, T, u, v))[0]
    return result

# Define the parameters
lbda = 10
T = np.pi
t = np.linspace(0,T,N)

n=0
while n<N:
  u = [integral_u(a, lbda, T, u0[n], v0[n]) for a in t]
  v = [integral_v(a, lbda, T, u0[n], v0[n]) for a in t]
  u0 = u
  v0 = v
  n = n+1

plt.plot(t,u)
plt.plot(t,v)
plt.show()

Answer 1

What you are doing is essentially a fixed-point iteration. That is, you are computing $\lim_{n\to\infty} U_n$, where $U_{n+1} = F(U_n)$. Here $U_n = (u_n,v_n)^T$ and $F$ is a vector containing your integral functions. Some typical convergence criteria for these types of problems is a combination of absolute and relative tolerances on the change in $U$ along with some maximum number of allowed iterations. For example, terminating the loop if any of the following conditions are met $$ \|U_{n+1} - U_n\|<\tau_a,\quad\frac{\|U_{n+1} - U_{n}\|}{\|U_1-U_0\|}<\tau_r,\quad n>N_{max} $$ The choice of norm here is important for function-valued iterates. For the infinite-dimensional setting, one typically takes something like $$\|U\|^2 = \int_0^T (u(t))^2dt + \int_0^T (v(t))^2dt.$$ Since you have discretized $u$ and $v$, you have the option of either translating these integrals into quadrature or computing the standard Euclidean norm on the discretized vectors themselves. Just to be explicit, you have the choice between $$\|U_n^h\|_2^2 = \sum_{i=1}^N u_{n,i}^2 + v_{n,i}^2 \ \text{or} \ \|U_n^h\|_q^2 = \sum_{i=1}^N w_i(u_{n,i}^2 + v_{n,i}^2),$$ where $w_i$ are some quadrature weights, e.g., trapezoid rule. Here the superscript $\cdot^h$ denotes that the function is discretized but I also use it to denote the uniform step size. For a problem this small the difference in convergence behavior is likely negligible, but one can check how they compare against the infinite-dimensional norm. $\|U_n^h\|_q$ will match the true norm to whatever the order of the quadrature scheme is, and $h\cdot\|U_n^h\|_2$ will match the true norm up to order $\mathcal{O}(h)$ for smooth functions $u$ and $v$.

Answer 2

You are computing convolutions. These are best evaluated as convolutions, rather than as sequences of integrals. A simple way of doing this is to discretize and then end up with a matrix-vector product where the matrix corresponds to the convolution operation.