I'm trying to solve the problem

enter image description here

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 :
    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)

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



2 Answers 2


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$.


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.

  • $\begingroup$ Even if the kernel is of the form $K(x-y)$ so it is a true convolution, the functions $u$ and $v$ are still inside some nonlinearity in the integral so some sort of nonlinear solver is needed $\endgroup$
    – whpowell96
    Commented May 18, 2023 at 17:21
  • $\begingroup$ I do agree that the step of performing the convolution is better handled as a matvec but it will need to go into some nonlinear solver in the end $\endgroup$
    – whpowell96
    Commented May 18, 2023 at 17:36
  • $\begingroup$ Sure, if you want to compute the limit you either need to iterate things out, or compute the solution of a nonlinear integral equation of the second kind. $\endgroup$ Commented May 18, 2023 at 23:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.