Fredholm第二类积分方程的Python代码实现（3）

这是《The Numerical Solution of Integral Equations of the Second Kind》书上的例题代码实现，P59-61。如有您发现有问题，欢迎指出，，，

# -*- coding: utf-8 -*-
"""
Created on Tue May 29 12:37:44 2018

@author: shaowu

本代码主要运用collocation 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)cos(t)和sqrt(t)，求出其对应的y(t)，然后再来反解x(t).更详细说明可
参见《The Numerical Solution of Integral Equations of the Second Kind》P59-61.
"""
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 gauss_xw(m=10):
    """
    默认用10个点求高斯——勒让德节点xi和权weight，并返回x和w数组
    """
    x,w=p_roots(m+1)
    return x,w

def gauss_solve_y(x,w,lamda,a,b,n): #参数n>=1
    """
    求解课本3.2.41式中的右端项y(ti).
    method:gauss_Legendre
    参数：
    x,w为高斯点和对应的权,
    a,b 对应于区间[a,b],
    n 子区间的个数
    返回的是一个n+1维的列表
    """
    c=(b-a)/2
    s=(b-a)/2*x+(a+b)/2 #把区间[a,b]变化到[-1,1]
    h=(b-a)/(n)
    t=[a+i*h for i in range(0,n+1)] #等距划分a=t0<t1<...<tn=b
    return [(lamda*scp.exp(-i)*scp.cos(i)-
            sum([c*w[j]*scp.exp(s[j]*i)*scp.exp(-s[j])*scp.cos(s[j]) for j in range(len(s))]))
            for i in t],\
           [(lamda*scp.sqrt(i)-
            sum([c*w[j]*scp.exp(s[j]*i)*scp.sqrt(s[j]) for j in range(len(s))]))
            for i in t]

def hat_f(j,x,a,b):
    """
    定义hat函数
    """
    h=(b-a)/(n)
    t=[a+i*h for i in range(0,n+1)]
    if j==0:
        if x>=t[0] and x<=t[1]:
            return 1-abs(x-t[0])/h
        else:
            return 0
    if j==n:
        if t[n-1]<=x and x<=t[n]:
            return 1-abs(x-t[n])/h
        else:
            return 0
    if j>0 and j<n:
        if t[j-1]<=x and x<=t[j+1]:
            return 1-abs(x-t[j])/h
        else:
            return 0
def elements_of_matrix(x,w,i,j,t,a,b,n):
    """
    对K(t,s)*lj(s)在[a,b]上求积分,即等式（3.2.41）中的积分项；
    返回的是一个数，即对应系数矩阵中的一个元素
    """
    h=(b-a)/(n)
    if j==0:
        #c=h/2;s=h/2*x+(t[0]+t[1])/2 #区间变化
        return -scp.integrate.quad(lambda x: (x-t[0])*scp.exp(t[i]*x),t[0],t[1])[0]/h
        #return -sum([w[k]*c*hat_f(j,s[k],a,b)*scp.exp(t[i]*s[k]) for k in range(len(s))])
    elif j==n:
        #c=h/2;s=h/2*x+(t[n-1]+t[n])/2 #区间变化
        return -scp.integrate.quad(lambda x: (t[n]-x)*scp.exp(t[i]*x),t[n-1],t[n])[0]/h
        #return -sum([w[k]*c*hat_f(j,s[k],a,b)*scp.exp(t[i]*s[k]) for k in range(len(s))])
    else:
        #c=h;s=h*x+(t[j-1]+t[j+1])/2 #区间变化
        return -scp.integrate.quad(lambda x: (x-t[j-1])*scp.exp(t[i]*x),t[j-1],t[j])[0]/h-\
                scp.integrate.quad(lambda x: (t[j+1]-x)*scp.exp(t[i]*x),t[j],t[j+1])[0]/h
        #return -sum([w[k]*c*hat_f(j,s[k],a,b)*scp.exp(t[i]*s[k]) for k in range(len(s))])
def solve_xn(x,w,lamda,a,b,n,y1,y2):
    """
    计算xn1和xn2
    返回xn1和xn2(一个n+1维的列表)
    """
    xn=[]
    h=(b-a)/(n)
    t=[a+i*h for i in range(0,n+1)] #等距划分a=t0<t1<...<tn<tn+1=b
    #A.append([-scp.integrate.quad(lambda x:(1-abs(x-t[j])/h)*scp.exp(t[0]*x),t[0],t[1])[0] for j in range(n+1)]) #计算i=0时的
    for i in range(0,n+1):
        #guodu=[-scp.integrate.quad(lambda x:(1-abs(x-t[0])/h)*scp.exp(t[i]*x),t[0],t[1])[0]]+\
        guodu=[elements_of_matrix(x,w,i,j,t,a,b,n) for j in range(0,n+1)]#+\
               #   [-scp.integrate.quad(lambda x:(1-abs(x-t[n])/h)*scp.exp(t[i]*x),t[n-1],t[n])[0]]
        xn.append(guodu)
    #A.append([-scp.integrate.quad(lambda x:(1-abs(x-t[j])/h)*scp.exp(t[n]*x),t[n-1],t[n])[0] for j in range(n+1)])  #计算i=n时的
    xn=np.array(xn)

    for i in range(n+1):
        xn[i][i]=lamda+xn[i][i]
    
    result1=np.linalg.solve(xn,y1) #求解xn1
    result2=np.linalg.solve(xn,y2) #求解xn2
    return result1,result2
def solve_error(lamda,a,b,n,xn1,xn2):
    """
    计算最大误差，并返回
    """
    h=(b-a)/(n)
    t=np.array([a+i*h for i in range(0,n+1)])
    #tt=np.array([i for i in np.arange(a+0.1,b,0.1)])
    xt1=scp.exp(-t)*scp.cos(t)  ##第一个精确解xt1
    xt2=scp.sqrt(t)   ##第二个精确解xt2
    #E1=LA.norm(xt1-xn,np.inf)  ##计算无穷范数
    E1=max(abs(xt1-xn1))
    #E2=LA.norm(xt2-xnt2,np.inf)
    E2=max(abs(xt2-xn2))
    print('The maximum error on the %d\'s collocation node points for xt1 is %s'%(n,E1))
    print('The maximum error on the %d\'s collocation node points for xt2 is %s'%(n,E2))
    """
    #画图
    plt.figure(1)
    plt.subplot(211)
    plt.title('Figure1.(where a=0,b=1,n=128,lambda=5)')
    plt.xlabel('t')
    plt.ylabel('x')
    plt.plot(t,xt1,'*--',t,xn1,'-')
    plt.subplot(212)
    plt.plot(t,xt2,'--+',t,xn2,'-')
    """
    return E1,E2
if __name__ == '__main__':
    print('******************************程序入口*******************************')
    lamda=int(input('pleas input lambda:'))
    n=int(input('please input n:')) # n>0 为配置点个数
    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)
    y1,y2=gauss_solve_y(x,w,lamda,a,b,n)
    xn1,xn2=solve_xn(x,w,lamda,a,b,n,y1,y2)
    E1,E2=solve_error(lamda,a,b,n,xn1,xn2)
    print('all_cost_time:',time.time()-start_time)

运行结果：