Fredholm第二类积分方程的Python代码实现-乘积积分法

 本代码主要运用Product integration method实现Fredholm integral equations of the second kind的求解,
参见《The Numerical Solution of Integral Equations of the Second Kind》P116-118.

# -*- coding: utf-8 -*-
"""
Created on Thu May 31 15:37:03 2018

@author: shaowu


本代码主要运用Product integration method实现Fredholm integral equations of the second kind的求解,方程如下：
    lamda*x(t)-integrate(K(t,s)*x(s))ds=y(t),a=0<=t<=b,K(t,s)取exp(s*t)
首先，我们给定精确解exp(t)，求出其对应的y(t)，然后再来反解x(t).更详细说明可
参见《The Numerical Solution of Integral Equations of the Second Kind》P116-118.
"""
import sympy as sp
import scipy as scp
import numpy as np
import pandas as pd
import numpy.linalg as LA
import matplotlib.pyplot as plt
from scipy.special.orthogonal import p_roots
import time
start_time=time.time()

def function_x(x):
    """
    定义精确解x
    """
    return scp.exp(x)

def function_k(t,s):
    """
    定义核函数K
    """
    return scp.log(abs(s-t))

def function_L(t,s):
    """
    定义L(t,s)函数
    """
    return 1

def gauss_xw(m=400):
    """
    默认用400个点求高斯——勒让德节点xi和权weight，并返回x和w都是一个数组
    """
    x,w=p_roots(m+1)
    return x,w

def gauss_solve_y(x,w,lamda,a,b,n):#参数x,w为高斯点和对应的权
    """
    用Gauss_Legendre积分公式求解方程的右端项y
    """
    h=(b-a)/(n)
    t=[a+i*h for i in range(0,n+1)] #等距划分a=t0<t1<...<tn=b
    y=[]
    for i in t:
        c1=(i-a)/2
        s1=(i-a)/2*x+(a+i)/2 ##对区间做变换：[a,i]-->[-1,1]
        c2=(b-i)/2
        s2=(b-i)/2*x+(i+b)/2 ##对区间做变换：[i,b]-->[-1,1]
        if i==a:
            y.append(lamda*function_x(i)-\
            sum([c2*w[j]*scp.log(s2[j]-i)*function_x(s2[j]) for j in range(len(s2))]))
        elif i==b:
            y.append(lamda*function_x(i)-
            sum([c1*w[j]*scp.log(i-s1[j])*function_x(s1[j]) for j in range(len(s1))]))
        else:
            y.append(lamda*function_x(i)-
            sum([c1*w[j]*scp.log(i-s1[j])*function_x(s1[j]) for j in range(len(s1))])-\
            sum([c2*w[j]*scp.log(s2[j]-i)*function_x(s2[j]) for j in range(len(s2))]))
    return y

def gauss_int1(x,w,i,lamda,a,b):#这是积分上限减去自变量的情况，即被积函数形式为（b-x）*f(x)
    """
    用Gauss_Legendre积分公式求权函数w中的积分项，课本(4.2.66)式中的积分项,被积函数形式为（b-x）*f(x)
    """
    c=(b-a)/2
    s=(b-a)/2*x+(a+b)/2 ##对区间做变换：[a,i]-->[-1,1]
    if i<=a:
        return sum([c*w[j]*(b-s[j])*scp.log(s[j]-i) for j in range(len(s))])
    else:
        return sum([c*w[j]*(b-s[j])*scp.log(i-s[j]) for j in range(len(s))])

def gauss_int2(x,w,i,lamda,a,b):#这是自变量减去积分下限的情况，即被积函数形式为（x-a）*f(x)
    """
    用Gauss_Legendre积分公式求权函数w中的积分项,被积函数形式为（x-a）*f(x)
    """
    c=(b-a)/2
    s=(b-a)/2*x+(a+b)/2 ##对区间做变换：[a,i]-->[-1,1]
    if i<=a:
        return sum([c*w[j]*(s[j]-a)*scp.log(s[j]-i) for j in range(len(s))])
    else:
        return sum([c*w[j]*(s[j]-a)*scp.log(i-s[j]) for j in range(len(s))])

def solve_weight(x,w,i,lamda,a,b,n): 
    """
    这里是计算权函数w(i)
    返回w是一个list列表
    """
    h=(b-a)/(n)
    t=[a+k*h for k in range(0,n+1)] #等距划分a=t0<t1<...<tn=b
    weight=[]
    for j in range(n+1):
        if j==0:
            weight.append(gauss_int1(x,w,i,lamda,t[0],t[1])/h)
        elif j==n:
            weight.append(gauss_int2(x,w,i,lamda,t[n-1],t[n])/h)
        else:
            weight.append(gauss_int2(x,w,i,lamda,t[j-1],t[j])/h+\
                          gauss_int1(x,w,i,lamda,t[j],t[j+1])/h)
    return weight

def solve_xn(x,w,lamda,a,b,n,y):
    """
    求解系数矩阵，并求出xn
    """
    A=[] #用于存放系数矩阵
    h=(b-a)/(n)
    t=[a+i*h for i in range(0,n+1)] #等距划分a=t0<t1<...<tn<tn+1=b
    for i in t:
        A.append(solve_weight(x,w,i,lamda,a,b,n)) #调用solve_weight函数
    A=np.array(A)
    A=-A
    for i in range(0,n+1):
        A[i][i]=lamda+A[i][i]
    xn=np.linalg.solve(A,y)
    return xn

def solve_error(xn,a,b,n):
    """
    计算误差
    """
    h=(b-a)/(n)
    t=[a+i*h for i in range(0,n+1)] #等距划分a=t0<t1<...<tn<tn+1=b
    E=LA.norm(function_x(t)-xn,np.inf)
    return E

if __name__ == '__main__':
    print('******************************程序入口*******************************')
    lamda=int(input('pleas input lambda:'))
    n=int(input('please input  n:'))
    a=int(input('please input the left value of interval:'))
    b=int(input('please input the right value of interval:'))
    m=int(input('please input the node of Gauss_Legendre:'))
    print('计算中...')
    x,w=gauss_xw(m)
    y=gauss_solve_y(x,w,lamda,a,b,n)
    xn=solve_xn(x,w,lamda,a,b,n,y)
    E=solve_error(xn,a,b,n)
    print('the error is:',E)
    print('all_cost_time:',time.time()-start_time)