数值计算/数值分析 python实践

import matplotlib.pyplot as plt
import numpy as np
from pylab import mpl
import math
import random
from sympy import *

# 通用的基函数
def LagrangeInterpolationDotN_Li(x, x_data, y_data, k):
    size = len(x_data)
    i = 0
    
    Ly = y_data[k]  # 初值为Yk
    
    while( i < size):
        if(i != k):
            Ly = Ly * (x-x_data[i])/(x_data[k]-x_data[i])
        i += 1
    return (Ly)
# 通用的插值函数
def LagrangeInterpolationDotN_F(x, x_data, y_data):
    size = len(x_data)
    k = 0
    sum = 0  # 初值为0
 
    while(k < size):
        sum = sum + LagrangeInterpolationDotN_Li(x, x_data, y_data, k)
        k += 1
    return (sum)
# 第一题需要用到的fx= 1 / (1 + x^2） 这个公式
def f(x):
    return 1 / (1 + x ** 2)
#第二题需要用的被积函数e^(2x)
def f2(x):
    return math.exp(2*x)
'''
1.1 第一题第一问：
'''
if __name__=="__main__":
    print('令插值节点为等距节点：-5、-4、-3、-2、-1、0、1、2、3、4、5，在这些节点处对 f (x)进行拉格朗日插值；')
    x_samples_dot11 = [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5]#创建x列表存储数据x值
    y_samples_dot11 = [1/(1+i**2) for i in x_samples_dot11] #创建y列表存储数据的y值 这里使用了列表推导式！
    print("X:", x_samples_dot11)
    print("Y:", y_samples_dot11)

    x_data_predict_dot11 = np.arange(-5,5.1,0.1) # 预测点的位置
    y_data_predict_dot11_F  = LagrangeInterpolationDotN_F (x_data_predict_dot11, x_samples_dot11, y_samples_dot11)
    y_data_real = [1/(1+i**2) for i in x_data_predict_dot11]
    error = np.abs(y_data_real - y_data_predict_dot11_F)
    max_error = np.max(error)
    print('误差值是：', max_error)
    plt.scatter(x_samples_dot11, y_samples_dot11, label="sample", color="black")
    plt.plot   (x_data_predict_dot11, y_data_predict_dot11_F, label="f(x)", color="pink")
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    plt.title("第一题第一问-拉格朗日插值拟合过程")
    plt.legend(loc="upper left")
    plt.show()

'''
1.2 第一题第二问：
'''
if __name__=="__main__":
    print('令插值节点为区间[-5,5]上的 11 次切比雪夫多项式的零点，在这些节点处对 f (x)进行拉格朗日插值；')
    number = [1,2,3,4,5,6,7,8,9,10,11]
    x_samples_Chebyshev = [0 + 5*math.cos((2*i - 1)*math.pi/(2*11)) for i in number]#创建x列表存储数据x值
    y_samples_Chebyshev = [1/(1+i**2) for i in x_samples_Chebyshev] #创建y列表存储数据的y值 这里使用了列表推导式！
    print("X:", x_samples_Chebyshev)
    print("Y:", y_samples_Chebyshev)

    x_data_predict_Chebyshev = np.arange(-5,5.1,0.1) # 预测点的位置
    y_data_predict_Chebyshev_F  = LagrangeInterpolationDotN_F (x_data_predict_Chebyshev, x_samples_Chebyshev, y_samples_Chebyshev)
    y_data_predict_Chebyshev_real = [1/(1+i**2) for i in x_data_predict_Chebyshev]
    print('误差值是：', np.max(np.abs(y_data_predict_Chebyshev_real - y_data_predict_Chebyshev_F)))
    plt.scatter(x_samples_Chebyshev, y_samples_Chebyshev, label="sample", color="black")
    plt.plot   (x_data_predict_Chebyshev, y_data_predict_Chebyshev_F, label="f(x)", color="orange")
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    plt.title("第一题第二问-拉格朗日插值拟合过程")
    plt.legend(loc="upper left")
    plt.show()

'''
1.3 第一题第三问：
'''
if __name__=="__main__":
    print('令插值节点为等距节点：-5、-4、-3、-2、-1、0、1、2、3、4、5，在这些节点处对 f (x)进行分段线性插值；')
    x_samples_1_3 = [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5]#创建x列表存储数据x值
    y_samples_1_3 = [1/(1+i**2) for i in x_samples_1_3]#创建y列表存储数据的y值
    print("X:", x_samples_1_3)
    print("Y:", y_samples_1_3)
    x_data_predict_1_3 = np.arange(-5,5.1,0.1) # 预测点的位置
    y_data_predict_1_3 = np.interp(x_data_predict_1_3 , x_samples_1_3, y_samples_1_3)
    y_data_predict_1_3_real = [1/(1+i**2) for i in x_data_predict_1_3]
    print('误差值是：', np.max(np.abs( y_data_predict_1_3_real - y_data_predict_1_3)))
    plt.scatter(x_samples_1_3, y_samples_1_3, label="sample", color="black")#画点
    plt.plot   (x_samples_1_3, y_samples_1_3, label="sample", color="pink")#画线
    plt.title("第一题第三问-分段线性插值")
    plt.legend(loc="upper left")
    plt.show()

'''
1.4 第一题第四问：
'''
print('令插值节点为等距节点：-5、-4、-3、-2、-1、0、1、2、3、4、5，在这些节点处对 f (x)进行三次样条插值：/n插值函数为 S(x) ，其中边界条件为第一类边界条件')

def cal(begin, end, i):
    by = f(begin)
    ey = f(end)
    I = Ms[i] * ((end - n) ** 3) / 6 + Ms[i + 1] * ((n - begin) ** 3) / 6 + (by - Ms[i] / 6) * (end - n) + (
            ey - Ms[i + 1] / 6) * (n - begin)
    return I

def ff(x):  # f[x0, x1, ..., xk]
    ans = 0
    for i in range(len(x)):
        temp = 1
        for j in range(len(x)):
            if i != j:
                temp *= (x[i] - x[j])
        ans += f(x[i]) / temp
    return ans

def calM():
    lam = [1] + [1 / 2] * 9
    miu = [1 / 2] * 9 + [1]
    # Y = 1 / (1 + n ** 2)
    # df = diff(Y, n)
    x = np.array(range(11)) - 5
    # ds = [6 * (ff(x[0:2]) - df.subs(n, x[0]))]
    ds = [6 * (ff(x[0:2]) - 1)]
    for i in range(9):
        ds.append(6 * ff(x[i: i + 3]))
    # ds.append(6 * (df.subs(n, x[10]) - ff(x[-2:])))
    ds.append(6 * (1 - ff(x[-2:])))
    Mat = np.eye(11, 11) * 2
    for i in range(11):
        if i == 0:
            Mat[i][1] = lam[i]
        elif i == 10:
            Mat[i][9] = miu[i - 1]
        else:
            Mat[i][i - 1] = miu[i - 1]
            Mat[i][i + 1] = lam[i]
    ds = np.mat(ds)
    Mat = np.mat(Mat)
    Ms = ds * Mat.I
    return Ms.tolist()[0]

def calnf(x):
    nf = []
    for i in range(len(x) - 1):
        nf.append(cal(x[i], x[i + 1], i))
    return nf

def calf(f, x):
    y = []
    for i in x:
        y.append(f.subs(n, i))
    return y

def nfSub(x, nf):
    tempx = np.array(range(11)) - 5
    dx = []
    for i in range(10):
        labelx = []
        for j in range(len(x)):
            if x[j] >= tempx[i] and x[j] < tempx[i + 1]:
                labelx.append(x[j])
            elif i == 9 and x[j] >= tempx[i] and x[j] <= tempx[i + 1]:
                labelx.append(x[j])
        dx = dx + calf(nf[i], labelx)
    return np.array(dx)

def draw(nf):
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    x = np.linspace(-5, 5, 101)
    y = f(x)
    Ly = nfSub(x, nf)
    plt.plot(x, y, label='原函数')
    plt.plot(x, Ly, label='三次样条插值函数')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()

    plt.savefig('1.png')
    plt.show()

def lossCal(nf): # 计算误差
    x = np.linspace(-5, 5, 101)
    y = f(x)
    Ly = nfSub(x, nf)
    Ly = np.array(Ly)
    temp = np.max(np.abs(Ly - y))  # 要求的误差
    print('误差值是：',temp)
    return temp
if __name__=="__main__":
    x = np.array(range(11)) - 5  # 根据题意-5~5
    y = f(x)
    n, m = symbols('n m')
    init_printing(use_unicode=True)
    Ms = calM()
    nf = calnf(x)
    draw(nf)
    lossCal(nf)

'''
1.5 第一题第五问：
'''
print('考虑等距节点：-5、-4、-3、-2、-1、0、1、2、3、4、5，求 f (x)在这些节点上的六次最小二乘拟合多项式,权重均为 1,')

class LeastSqauresCF():
    def __init__(self):
        self.order = 1
        self.Xvalues = None
        self.Yvalues = None
        self.fig = plt.figure()
        self.subfig = self.fig.add_subplot(111)
    #========================================
    # Draw the input points/curve in a figure
    #========================================
    def ReadAnddrawInputPoints(self,Inputx=None,InputY=None, order=1):
        try:
            if Inputx is None or InputY is None:
                return False
            if(len(Inputx)!=len(InputY)):
                return False
            if order <=0:
                return False
            
            if type(order) !=int:
                self.order = int(order)
            else:
                self.order = order
            self.Xvalues = Inputx
            self.Yvalues = InputY
            
            self.subfig.plot(Inputx,InputY,color='m',linestyle='',marker='*')
            return True
        except Exception as e:
            print(e)
            return False
    #========================================
    # Calculate curve fitting factor matrix
    #========================================    
    def produceFittingCurveFactors(self):
        try:
            #(1) Generate matrix [X]
            matX=[]
            for i in range(0,len(self.Xvalues)):
                matx1=[]
                for j in range(0,self.order+1):
                    dx=1.0
                    for l in range(0,j):
                        dx = dx * self.Xvalues[i]
                    matx1.append(dx)
                matX.append(matx1)        

            #(2) Generate matrix [X]T.[X]
            matX_Trans = np.matrix(matX).T
            matX_FinalX = np.dot(np.matrix(matX_Trans),np.matrix(matX))

            #(3) Generate matrix Y' =[X]T.[Y]
            matFinalY = np.dot(matX_Trans,np.matrix(self.Yvalues).T)

            #(4) Solve the function:[A] = [[X]T.[X]]**(-1).[X]T.[Y]
            matAResult=np.linalg.solve(np.array(matX_FinalX),np.array(matFinalY))
            return matAResult
        except Exception as e:
            print(e)
            return None
    #========================================
    # Output fitting curve function
    #========================================        
    def outputFittingCurveFunction(self, inputMatFactors=None):
        if inputMatFactors is None:
            return False
        i = 0
        strFitting="Function(x)="
        for a in inputMatFactors:
            #print(a[0])
            if i==0:
                strFitting +=str(a[0])
            else:
                strFitting +=("+"if a[0]>0 else"") +str(a[0])+(("x**"+str(i)) if i>1 else "x")
            i+=1
        print(strFitting)
        return strFitting
    #========================================
    # draw the curve based on the result function
    #========================================
    def drawFittedCurve(self, xRangeMin=None,xRangeMax=None,matAResult=None,resolution=0.01):
        try:
            if matAResult is None:
                return False
            if xRangeMin is None or xRangeMax is None:
                xRangeMin = self.Xvalues[0]
                xRangeMax = self.Xvalues[-1]
            
            #print('xRangeMin: ',xRangeMin, 'xRangeMax: ',xRangeMax)
            
            xxa= np.arange(xRangeMin,xRangeMax,resolution)
            yya=[]
            for i in range(0,len(xxa)):
                yy=0.0
                for j in range(0,self.order+1):
                    dy=1.0
                    #x[i]**j
                    for k in range(0,j):
                        dy*=xxa[i]
                    #a[j]*(x[i]**j)
                    dy*=matAResult[j]
                    yy+=dy
                yya.append(yy)
            #print(xxa,yya)
            self.subfig.plot(xxa,yya,color='g',linestyle='-',marker='')
            plt.show() 
            return True
        except Exception as e:
            print(e)
            return False
    def calculateAbsoluteError(self, matAResult):
        try:
            yFitted = []
            for x in self.Xvalues:
                y = 0.0
                for i, a in enumerate(matAResult):
                    y += a * (x ** i)
                yFitted.append(y)
            absoluteError = np.abs(np.array(self.Yvalues) - np.array(yFitted))
            meanAbsoluteError = np.mean(absoluteError)
            return meanAbsoluteError
        except Exception as e:
            print(e)
            return None

if __name__=="__main__":
    LS = LeastSqauresCF()
    sorted_x = [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5]
    sorted_y = [1/(1+i**2) for i in x_samples_dot11]
    if LS.ReadAnddrawInputPoints(sorted_x,sorted_y,6):	
        MatrixFactor = LS.produceFittingCurveFactors()
        LS.outputFittingCurveFunction(MatrixFactor)
        LS.drawFittedCurve(matAResult=MatrixFactor)
        print("误差值是:",LS.calculateAbsoluteError(MatrixFactor))

'''
2.1 第二题第一问：
'''
print('将区间[3,5]等分 40 份，用复合梯形公式计算:')
def trapezoid(a, b, n): # 复合梯形公式
    h=(b-a)/n
    x=a
    s=f2(x)-f2(b)
    for k in range(1,n+1):
        x=x+h
        s=s+2*f2(x)
    result=(h/2)*s
    return result

def simpson(a,b,n):   # 辛普森公式
    h=(b-a)/n
    x=a
    s=f2(x)-f2(b)
    for k in range(1,n+1):
        x=x+h/2
        s=s+4*f2(x)
        x=x+h/2
        s=s+2*f2(x)
    result=(h/6)*s
    return result
    
if __name__=="__main__":
    print('S=',trapezoid(3,5,40))

'''
2.2 第二题第二问：
'''
print('将区间[3,5]等分 20 份，用复合辛普森公式计算（每个小区间需加中点）：')
if __name__=="__main__":
    print('S=',simpson(3,5,40))

'''
2.3 第二题第三问：
'''
print('T0将区间[3,5]等分为 20 份的复合梯形公式，用龙贝格积分算法加速一次,算出T1：')
import math


def romberg_integration(a, b, n, num_iterations,f):
    results = [[0] * (num_iterations+1) for _ in range(num_iterations+1)]
    h = (b - a) / n
    x = [a + i * h for i in range(n+1)]
    sum_value = sum([f(x[i]) for i in range(1, n)])
    results[0][0] = h * (0.5 * (f(a) + f(b)) + sum_value)
    
    for i in range(1, num_iterations+1):
        h *= 0.5
        sum_value = sum([f(x[j]) for j in range(1, 2**(i-1)+1, 2)])
        results[i][0] = 0.5 * results[i-1][0] + h * sum_value
        
        for j in range(1, i+1):
            results[i][j] = results[i][j-1] + (results[i][j-1] - results[i-1][j-1]) / ((4**j) - 1)
    
    return results[num_iterations][num_iterations]

if __name__=="__main__":
    # 使用龙贝格积分算法加速一次
    T1_h_romberg = romberg_integration(3, 5, 20, 1, f2)
    print("T1(h) (with Romberg acceleration) =", T1_h_romberg)

'''
2.4 第二题第四问：
'''
print('将区间[3,5]等分为 10 份的复合梯形公式，用龙贝格积分算法加速两次,算出T2：')
if __name__=="__main__":
    T2_h_romberg = romberg_integration(3, 5, 10, 2, f2)
    print("T2(h) (with Romberg acceleration) =", T2_h_romberg)

'''
2.5 第二题第五问：
'''
print('将区间[3,5]等分为 5 份的复合梯形公式，用龙贝格积分算法加速三次,算出32：')
if __name__=="__main__":
    T3_h_romberg = romberg_integration(3, 5, 5, 3, f2)
    print("T3(h) (with Romberg acceleration) =", T3_h_romberg)

'''
2.6 第二题第六问：
'''
print('将区间[3,5]等分 20 份，用复合的 2 点高斯公式计算：')

def composite_gauss2(f, a, b, n):
    # 复合的2点高斯公式计算积分 书上有结论~
    integral = 0
    nodes = np.linspace(a, b, n+1)
    
    for i in range(n):
        x1 = nodes[i]
        x2 = nodes[i+1]
        mid_point = 0.5 * (x1 + x2)
        width = 0.5 * (x2 - x1)
        
        # 计算节点上的权重和函数值
        weight1 = 1 # A0=A1=1  X0 X1是正负跟号3
        weight2 = 1
        function_value1 = f(0.5 * width * (-1/np.sqrt(3)) + mid_point) # 平移到子区间
        function_value2 = f(0.5 * width * (1/np.sqrt(3)) + mid_point)
        
        # 计算积分部分并累加
        integral += width * (weight1 * function_value1 + weight2 * function_value2)
    
    return integral

# 将区间[3, 5]等分为20份，计算积分
integral01 = composite_gauss2(f2, 3, 5, 20)
print("Integral:", integral01)

'''
2.7 第二题第七问：
'''
print('将区间[3,5]等分 10 份，用复合的 4 点高斯公式计算：')
def composite_gauss4(f, a, b, n):
    # 复合的四点高斯公式计算积分
    integral = 0
    nodes = np.linspace(a, b, n+1)
    
    for i in range(n):
        x1 = nodes[i]
        x2 = nodes[i+1]
        h = x2 - x1
        
        # 四个节点的位置
        t1 = -np.sqrt((3 + 2*np.sqrt(6/5))/7)
        t2 = -np.sqrt((3 - 2*np.sqrt(6/5))/7)
        t3 = np.sqrt((3 - 2*np.sqrt(6/5))/7)
        t4 = np.sqrt((3 + 2*np.sqrt(6/5))/7)
        
        # 四个节点的权重
        w1 = (18 - np.sqrt(30))/36
        w2 = (18 + np.sqrt(30))/36
        w3 = (18 + np.sqrt(30))/36
        w4 = (18 - np.sqrt(30))/36
        
        # 计算积分部分并累加
        integral += 0.5 * h * (w1 * f(0.5 * h * t1 + 0.5 * (x1 + x2)) +
                              w2 * f(0.5 * h * t2 + 0.5 * (x1 + x2)) +
                              w3 * f(0.5 * h * t3 + 0.5 * (x1 + x2)) +
                              w4 * f(0.5 * h * t4 + 0.5 * (x1 + x2)))
    
    return integral

# 将区间[3, 5]等分为20份，计算积分
integral02 = composite_gauss4(f2, 3, 5, 10)
print("Integral:", integral02)


'''
3 第三题：
'''
A = np.array([[3, 1, 1], [1, 3, 1], [1, 1, 3]])
b = np.array([3, 0, 2])
x0 = np.zeros_like(b)  # 初始猜测解 全0
max_iterations = 100   # 最大迭代次数
tolerance = 1e-8       # 允许误差

# 使用numpy.linalg.solve()函数求解线性方程组
solution_true = np.linalg.solve(A, b)
print("正确的解为", solution_true)

print('雅可比迭代法：')
def DLU(A):
    D=np.zeros(np.shape(A))
    L=np.zeros(np.shape(A))
    U=np.zeros(np.shape(A))
    for i in range(A.shape[0]):
        D[i,i]=A[i,i]
        for j in range(i):
            L[i,j]=-A[i,j]
        for k in list(range(i+1,A.shape[1])):
            U[i,k]=-A[i,k]
    L=np.mat(L)
    D=np.mat(D)
    U=np.mat(U)
    return D,L,U
 
#迭代
def Jacobi_iterative(A,b,x0,maxN,p):  #x0为初始值，maxN为最大迭代次数，p为允许误差
    D,L,U=DLU(A)
    if len(A)==len(b):
        D_inv=np.linalg.inv(D)
        D_inv=np.mat(D_inv)
        B=D_inv * (L+U)
        B=np.mat(B)
        f=D_inv*b
        f=np.mat(f)
    else:
        print('维数不一致')
        sys.exit(0)  # 强制退出
    
    a,b=np.linalg.eig(B) #a为特征值集合，b为特征值向量
    c=np.max(np.abs(a)) #返回谱半径
    if c<1:
        print('迭代收敛')
    else:
        print('迭代不收敛')
        sys.exit(0)  # 强制退出
#以上都是判断迭代能否进行的过程，下面正式迭代
    k=0
    while k<maxN:
        x=B*x0+f
        k=k+1
        eps=np.linalg.norm(x-x0,ord=2)
        if eps<p:
            break
        else:
            x0=x
    return k,x

A1 = np.mat([[3, 1, 1], [1, 3, 1], [1, 1, 3]])
b1 = np.mat([[3],[0],[2]])
x01=np.mat([[0],[0],[0]])
k,solution_jacobi=Jacobi_iterative(A1,b1,x01,max_iterations,tolerance)
print("迭代次数:",k)
residual = A1.dot(solution_jacobi) - b1
residual_norm = np.linalg.norm(residual, ord=np.inf)
print("残差：", residual_norm)
print("Solution:", solution_jacobi)

print('高斯赛德尔迭代法')
def Guss_Selbi(a, b, x, g):  # a为系数矩阵  b增广的一列  x迭代初始值  g计算精度
    x = x.astype(float)  # 设置x的精度，让x计算中能显示多位小数
    m, n = a.shape
    times = 0  # 迭代次数
    if (m < n):
        print("There is a 解空间。")  # 保证方程个数大于未知数个数
    else:
        while True:
            for i in range(n):
                s1 = 0
                tempx = x.copy()  # 记录上一次的迭代答案
                for j in range(n):
                    if i != j:
                        s1 += x[j] * a[i][j]
                x[i] = (b[i] - s1) / a[i][i]
                times += 1    # 迭代次数加一
            gap = max(abs(x - tempx))  # 与上一次答案模的差
 
            if gap < g:  # 精度满足要求，结束
                break
 
            elif times > 10000:  # 如果迭代超过10000次，结束
                break
                print("10000次迭代仍不收敛")
 
    print('迭代次数：',times)
    return x

solution_gauss = Guss_Selbi(A, b, x0, tolerance)
print("残差：", np.linalg.norm(solution_true - solution_gauss, ord=np.inf))
print("Solution:", solution_gauss)

print('超松弛迭代法：')
def sor_iteration(A,b,X0,max1, tolerance, w):
    '''A代表线性方程组AX=b的系数矩阵
           b代表线性方程组AX=b右边的部分
           X0代表高斯—赛德尔迭代的初始值
           w代表松弛因子
           max1代表迭代的次数'''
    n=np.shape(A)[0]
    L=-np.tril(A,-1)
    U=-np.triu(A,1)
    D=A+L+U
    B=np.dot(np.linalg.inv(D-w*L),((1-w)*D+w*U))
    f=np.dot(np.linalg.inv(D-w*L),w*b)
    err = np.inf
    for i in range(max1) :
        if err <= tolerance:
            break
        X0=np.dot(B,X0)+f
        err=np.linalg.norm(np.dot(A,X0)-b,ord=2)
    print("迭代次数：", i)
    return X0

omega = 6 / (3 + math.sqrt(5))
solution_sor = sor_iteration(A, b, x0, max_iterations, tolerance, omega)
print("残差：", np.linalg.norm(solution_true - solution_sor, ord=np.inf))
print("Solution:", solution_sor)

print('最速梯度下降法：')
def steepest_descent_linear_system(A, b, initial_x, max_iterations, tolerance):
    x = initial_x
    iteration = 0
    error = np.inf
    
    while iteration < max_iterations and error > tolerance:
        Ax = np.dot(A, x)
        residual = Ax - b
        gradient = np.dot(A.T, residual)
        alpha = np.dot(residual, residual) / np.dot(np.dot(A, residual), residual)
        
        x_new = x - alpha * residual
        error = np.linalg.norm(x_new - x)
        x = x_new
        iteration += 1
    print('迭代次数：', iteration)
    return x

solution_sdls = steepest_descent_linear_system(A, b, x0, max_iterations, tolerance)
print("残差：", np.linalg.norm(solution_true - solution_sdls, ord=np.inf))
print("Solution:", solution_sdls)