如何用Scipy实现积分方程组序列的收敛性编程

我需要求解如下迭代格式的积分方程：

$u_{n+1}(t) = \int_0^T K_\lambda(t,s)\left[\sigma(s) + \lambda u_n(s) + f(s,u_n(s),v_n(s))\right]ds$
$v_{n+1}(t) = \int_0^T K_\lambda(t,s)\left[\sigma(s) + \lambda v_n(s) + g(s,u_n(s),v_n(s))\right]ds$

其中$u_n$和$v_n$是收敛到解$u$和$v$的序列，$\lambda$、$\sigma$、$f$、$g$和$K_\lambda$均为已知量。

我尝试用Scipy的integrate模块结合while循环实现，但不清楚如何设置收敛条件，也不知道如何将$u_n(t)$作为被积函数。以下是我的尝试代码：

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

问题修正与优化方案 #

核心问题分析 #

• 原代码迭代时错误调用了u0[n]和v0[n]，应该传入整个序列而非当前迭代次数对应的单个值
• 缺少收敛判断逻辑，固定迭代次数无法保证解的收敛性
• 未处理离散序列在积分变量任意点的取值问题，无法将$u_n(t)$作为被积函数的一部分

修改后的代码 #

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

# 参数配置
N = 30
bta = 2.1
eps = 0.1
sgm = 1.5
lbda = 10
T = np.pi
t_nodes = np.linspace(0, T, N)

# 初始迭代序列
u_prev = np.ones(N) * 0.2
v_prev = np.ones(N) * 0.15

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

# 已知函数f和g
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

# 计算下一轮迭代的u序列
def compute_u_next(t_eval, u_prev, v_prev, lbda, T, t_nodes):
    u_next = np.zeros_like(t_eval)
    for idx, t_val in enumerate(t_eval):
        def integrand(s):
            # 用线性插值获取离散序列在积分点s处的值
            u_s = np.interp(s, t_nodes, u_prev)
            v_s = np.interp(s, t_nodes, v_prev)
            return K(t_val, s, lbda, T) * (sigma(s) + lbda * u_s + f(s, u_s, v_s))
        u_next[idx] = integrate.quad(integrand, 0, T)[0]
    return u_next

# 计算下一轮迭代的v序列
def compute_v_next(t_eval, u_prev, v_prev, lbda, T, t_nodes):
    v_next = np.zeros_like(t_eval)
    for idx, t_val in enumerate(t_eval):
        def integrand(s):
            u_s = np.interp(s, t_nodes, u_prev)
            v_s = np.interp(s, t_nodes, v_prev)
            return K(t_val, s, lbda, T) * (sigma(s) + lbda * v_s + g(s, u_s, v_s))
        v_next[idx] = integrate.quad(integrand, 0, T)[0]
    return v_next

# 迭代收敛控制
max_iter = 100
tol = 1e-6  # 收敛阈值
iter_count = 0
converged = False

while iter_count < max_iter and not converged:
    u_next = compute_u_next(t_nodes, u_prev, v_prev, lbda, T, t_nodes)
    v_next = compute_v_next(t_nodes, u_prev, v_prev, lbda, T, t_nodes)
    
    # 计算相邻迭代的最大误差
    error_u = np.max(np.abs(u_next - u_prev))
    error_v = np.max(np.abs(v_next - v_prev))
    
    if error_u < tol and error_v < tol:
        converged = True
        print(f"迭代收敛，迭代次数：{iter_count+1}")
    else:
        u_prev = u_next.copy()
        v_prev = v_next.copy()
        iter_count += 1

if not converged:
    print(f"达到最大迭代次数{max_iter}，未收敛")

# 结果可视化
plt.figure(figsize=(10,6))
plt.plot(t_nodes, u_prev, label='u(t)')
plt.plot(t_nodes, v_prev, label='v(t)')
plt.xlabel('t')
plt.ylabel('函数值')
plt.title('积分方程迭代解')
plt.legend()
plt.grid(True)
plt.show()

关键修改说明 #

• 被积函数处理：使用np.interp将离散的迭代序列插值到积分变量$s$的任意点，解决了离散序列作为被积函数的问题
• 收敛条件设置：通过计算相邻迭代序列的最大绝对误差，与预设阈值tol对比，满足条件则停止迭代
• 迭代逻辑修正：每次迭代传入完整的序列，保证积分时能获取任意$s$点的$u_n(s)$和$v_n(s)$值
• 稳定性优化：使用copy()避免数组引用导致的意外修改，确保迭代数据独立

内容的提问来源于stack exchange，提问作者tareqath