本部分运用Nystrom 方法求解Fredholm integral equations of the second kind,对应书上P100-103。
# -*- coding: utf-8 -*- """ Created on Wed May 30 21:16:58 2018 @author: shaowu 本代码主要运用Nystrom 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》P100-103. """ 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): """ 定义精确解 """ return scp.exp(x) def function_k(s,t): """ 定义核函数 """ return scp.exp(s*t) def simpson_solve_A(a,b,lamda,n): """ 用Simpson求积公式求解(4.1.5)式子中的系数矩阵 """ h=(b-a)/n t=[a+i*h for i in range(0,n+1)] A=[] for i in t: A.append([-(b-a)/6*function_k(a,i),-(b-a)/6*4*function_k((a+b)/2,i),\ -(b-a)/6*function_k(b,i)]) A=np.array(A) for i in range(A.shape[0]): A[i][i]=lamda+A[i][i] return A def simpsom_sovle_y(a,b,lamda,n): """ 由给定的精确解,计算Simpson点对应的右端项y """ h=(b-a)/n t=[a+i*h for i in range(0,n+1)] return [lamda*function_x(i)-scp.integrate.quad(lambda s:(function_k(s,i)*function_x(s)),a,b)[0]\ for i in t] def simpsom_solve_xn(A,y): """ 求解xn """ return np.linalg.solve(A,y) def simpsom_error(xn,a,b,n): """ 计算误差 """ h=(b-a)/n t=[a+i*h for i in range(0,n+1)] return function_x(t)-xn ######################################################################################################### def gauss_xw(n=2): """ 默认用3个点求高斯——勒让德节点xi和权weight,并返回x和w数组 """ x,w=p_roots(n+1) return x,w def gauss_sovle_y(x,a,b,lamda): """ 由给定的精确解,计算高斯点对应的右端项y """ t=(b-a)/2*x+(a+b)/2 #对区间[a,b]做变化 return [lamda*function_x(i)-scp.integrate.quad(lambda s:(function_k(s,i)*function_x(s)),a,b)[0]\ for i in t] def gauss_solve_A(x,w,a,b,lamda): """ 计算系数矩阵 """ c=(b-a)/2 t=(b-a)/2*x+(a+b)/2 #把区间[a,b]变化到[-1,1] A=[] for i in t: A.append([-c*w[j]*function_k(t[j],i) for j in range(len(t))]) A=np.array(A) for i in range(A.shape[0]): A[i][i]=lamda+A[i][i] return A def gauss_solve_xn(A,y): """ 计算xn """ return np.linalg.solve(A,y) def gauss_error(xn,x,a,b): """ 计算误差 """ t=(b-a)/2*x+(a+b)/2 return function_x(t)-xn if __name__ == '__main__': print('******************************程序入口*******************************') lamda=int(input('pleas input lambda:')) n=int(input('please input Simpson node 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('计算中...') #下面是使用3个Simpson节点(输入时n=2的情况): y=simpsom_sovle_y(a,b,lamda,n) A=simpson_solve_A(a,b,lamda,n) xn=simpsom_solve_xn(A,y) E1=simpsom_error(xn,a,b,n) print(E1) ############################################################################# #下面是使用m个高斯勒让德节点(输入时m相应少1): x,w=gauss_xw(m) y=gauss_sovle_y(x,a,b,lamda) A=gauss_solve_A(x,w,a,b,lamda) xn=gauss_solve_xn(A,y) E2=gauss_error(xn,x,a,b) print(E2) print('all_cost_time:',time.time()-start_time)
运行结果:
******************************程序入口*******************************
pleas input lambda:2
please input Simpson node n:2
please input the left value of interval:0
please input the right value of interval:1
please input the node of Gauss_Legendre:5
计算中...
[-0.0047043 -0.00796678 -0.01640422]
[ 2.72448730e-13 2.98427949e-13 3.51274565e-13 4.72955008e-13
8.05577827e-13 1.36868294e-12]
all_cost_time: 35.219531774520874