配置法 求解1D第二类线性的Fredholm积分方程

# -*- coding: utf-8 -*-
"""
    Created on Sat Oct 16 19:27:08 2021
    Collocation_1D_Second_Fredholm_IE
    
    这是一个采用配置法 求解第二类 线性的fredholm积分方程的程序
    
    基函数采用 分片的线性基函数/帽函数

    u(x)-Ku=f(x)
    
    Ku=k(x,y)*u(y)在[a,b]上关于y的积分

  
@author: Mr.Yang
"""
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la
def Gaussian_Integral(a,b,f):#定义一元函数的7点高斯积分区间在-1 到 1上
    x_y_w=np.array([[-0.9491079123427585,0.1294849661688697],
                    [-0.7415311855993945,0.2797053914892766],
                    [-0.4058451513773972,0.3818300505051189],
                    [ 0.0000000000000000,0.4179591836734694],
                    [ 0.4058451513773972,0.3818300505051189],
                    [ 0.7415311855993945,0.2797053914892766],
                    [ 0.9491079123427585,0.1294849661688697]])
    S=0
    for i in range(0,np.size(x_y_w,0)):
         xq = x_y_w[i,0]
         wq = x_y_w[i,1]
         S=S+wq*f(((b-a)/2)*xq+(b+a)/2)
    #print(S*(b-a)/2)
    return S*(b-a)/2
a1=0
b1=1
N=20
x=np.linspace(a1,b1,N)#将区间[a1,b1]剖分为 N个点

def K(x,y):#核函数
    f=x+y
    return f

def fn(x):#右端函数
    f=x/2-1/3
    return f

def e_fn(x):#精确解
    f=x
    return f

h=(b1-a1)/(N-1)#步长  均匀的

A=np.zeros((N,N))
F=np.zeros(N)#右端荷载向量  Ax=b中的b

I=np.eye(N,N)#单位阵

for i in range(0,N):
    F[i]=fn(x[i])
    for j in range(0,N):
        if(j==0):
            a=x[j]
            b=x[j+1]
            
            phi=lambda y:K(x[i],y)*(x[1]-y)/h
            A[i,j]=Gaussian_Integral(a,b,phi)
        elif(j==(N-1)):
            a=x[N-2]
            b=x[N-1]
            
            phi=lambda y:K(x[i],y)*(y-x[N-2])/h
            A[i,j]=Gaussian_Integral(a,b,phi)
        else:
            a=x[j-1]
            b=x[j]
            c=x[j+1]
            phil=lambda y:K(x[i],y)*(y-x[j-1])/h
            phir=lambda y:K(x[i],y)*(x[j+1]-y)/h
            A[i,j]=Gaussian_Integral(a,b,phil)+Gaussian_Integral(b,c,phir)
            #print(a,b,c)
#print(A.round(4))
#print(F)


u= la.solve(I-A,F)#计算得到没一点点的解
#print(u)
#%%
if e_fn is not None:
    #
    #  evaluate the exact solution at the nodes.
    #
    uex = np.zeros (N)
    err = np.zeros (N)
    for i in range (0, N):
        uex[i] = e_fn (x[i])
        err[i]= abs ( uex[i] - u[i] )
    #
    #  Compare the solution and the error at the nodes.
    #
    Show_Table={'Y':True, #显示误差表 Yes
                'N':False}#不显示误差表 No
    if Show_Table['Y']: #改变字母即可
        print ( '' )
        print ( '  Node          数值解           精确解          误差' )
        print ( '' )
        for i in range ( 0, N ):
            #err = abs ( uex[i] - u[i] )
            print ( '  %4d  %14.6g  %14.6g  %14.6g' % ( i, u[i], uex[i], err[i] ) )
    Show_Picture={'Y':True, #显示误差表 Yes
                  'N':False}#不显示误差表 No
    if(Show_Picture['Y']):
        npp = 51
        xp = np.linspace ( a1, b1, npp )
        up = np.zeros ( npp )
        for i in range ( 0, npp ):
            up[i] = e_fn ( xp[i] )
        
        plt.plot ( x, u, 'bo-', xp, up, 'r.' )
        #plt.savefig ( 'fem1d.png' )
        plt.show()
        
#%%
'''
Node          数值解           精确解          误差

     0    -9.69316e-17               0     9.69316e-17
     1       0.0526316       0.0526316     8.32667e-17
     2        0.105263        0.105263     9.71445e-17
     3        0.157895        0.157895     2.22045e-16
     4        0.210526        0.210526     1.11022e-16
     5        0.263158        0.263158     1.66533e-16
     6        0.315789        0.315789     1.11022e-16
     7        0.368421        0.368421     1.66533e-16
     8        0.421053        0.421053     3.33067e-16
     9        0.473684        0.473684     1.11022e-16
    10        0.526316        0.526316     1.11022e-16
    11        0.578947        0.578947     2.22045e-16
    12        0.631579        0.631579     3.33067e-16
    13        0.684211        0.684211     4.44089e-16
    14        0.736842        0.736842     1.11022e-16
    15        0.789474        0.789474     4.44089e-16
    16        0.842105        0.842105     3.33067e-16
    17        0.894737        0.894737               0
    18        0.947368        0.947368     2.22045e-16
    19               1               1     4.44089e-16
'''