本程式碼主要運用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)