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

So far this gives a result that still doesn't converge, but at least it looks correct and maybe the problem is ill-posed. If anyone knows any improvements feel free to add them

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

N = 50


u0 = np.ones(N)*0.10
v0 = np.ones(N)*0.12

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)

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)
err=1
epss=1e-6
n=0
while err>epss:
  u = [integral_u(a, lbda, T, np.interp(a,t,u0), np.interp(a,t,v0)) for a in t]
  v = [integral_v(a, lbda, T, np.interp(a,t,u0), np.interp(a,t,v0)) for a in t]  
  err=np.linalg.norm(np.array(u)-np.array(u0),2)+np.linalg.norm(np.array(v)-np.array(v0),2)
  err=np.sqrt(err) 
  u0 = u
  v0 = v
  n = n+1
  print('n = ',n)
  print('err = ',err)

  

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

print(u[0])