数值分析实验 第二章 插值方法

前言

期末了，要考试了，摸了一学期的鱼(听不懂)，现在都要偿还了。

主体

对于插值的理解(学渣的个人看法啊)

插值方法就是将复杂的函数f(x)简单化，找到f(x)的近似且为简单函数的函数p(x)

根据f(x)中已知的结点(xi,yi),i=0,1,…,n,可用多种方法找出其规律，并构造出近似且简单函数p(x)，一般情况下，p(x)的曲线与由已知结点连接而成的曲线越相似，越靠近，取不同的x0时，得到的结果p(x0)与真实值(x0,y0)之间的误差就越小

泰勒公式就是一种插值方法，可以说不学会它，数值分析就没法学明白(我就没学好)

Taylor公式（泰勒公式）通俗+本质详解

泰勒公式是做什么用的
用一个多项式函数去逼近一个给定的函数(即尽量使多项式函数图像拟合给定的函数图像)，注意，逼近的时候一定是从函数图像上的某个点展开。如果一个非常复杂函数，想求其某点的值，直接求无法实现，这时候可以使用泰勒公式去近似的求该值

1.拉格朗日插值多项式

借鉴文章
https://www.zhihu.com/question/58333118

import numpy as np

'''
定义了一个名为`lagrange_interpolation`的函数，它接受插值节点序列的x和y值以及要计算近似值的点x_interpolate。函数首先创建一个空列表L_values，用于存储每个Li(x)的值。然后，它遍历插值节点序列中的每个节点，计算Li(x)并将其添加到L_values列表中。

接下来，函数计算插值多项式的值，即每个Li(x)乘以相应的y值，然后将它们相加。最后，函数返回计算得到的近似值。

在示例中，我们使用了一个简单的插值节点序列，并计算了点0.5的近似值。
'''


def lagrange_interpolation(x_values, y_values, x_interpolate):
   n = len(x_values)
   L_values = []
   
   # 内循环 计算Li(x)
   for i in range(n):
       L = 1
       for j in range(n):
           if j != i:
               # li(x) -> (x-xj)(xk-xj)的从0到j(不为i)累积
               L *= (x_interpolate - x_values[j]) / (x_values[i] - x_values[j])
       L_values.append(L)
   
   # 外循环 求p(x)
   result = 0
   for i in range(n):
       #p(x) -> yk * li(x)
       result += y_values[i] * L_values[i]
   
   return result
'''
# 示例
x_values = [0, 1, 2, 3]
y_values = [1, 2, 3, 4]
x_interpolate = 0.5
'''

x_values = []
y_values = []

print("请输入插值节点的个数：")
n = int(input())

print("请输入插值节点的x值：")
for i in range(n):
   x = float(input())
   x_values.append(x)

print("请输入插值节点的y值：")
for i in range(n):
   y = float(input())
   y_values.append(y)

print("请输入要计算近似值的点x_interpolate：")
x_interpolate = float(input())

result = lagrange_interpolation(x_values, y_values, x_interpolate)
print("f(", x_interpolate, ") ≈ ", result)

成效

请输入插值节点的个数：
4
请输入插值节点的x值：
1
5
6
7
请输入插值节点的y值：
3
5
6
8
请输入要计算近似值的点x_interpolate：
3
f( 3.0 ) ≈  4.399999999999997

2.牛顿插值多项式

https://www.zhihu.com/question/22320408/answer/141973314

##注意：这里假设插值节点是按照x的顺序给出的。
def newton_interpolation(points, x):
    # 计算(x-x0)(x-x1)...(x-xi)
    def newton_basis(i):
        result = 1
        for j in range(i):
            result *= (x - points[j][0])
        return result

    #求差商
    def diff_quotient(i, j):
        if i == j-1:
            #f(xj)-f(xi)
            up=points[j][1]-points[i][1]
            #up/(xj-xi)
            return up/(points[i+1][0]-points[i][0])
        #递归
        up=diff_quotient(i+1,j)-diff_quotient(i,j-1)
        return up/(points[j][0]-points[i][0])
        
        
    # 计算牛顿插值多项式的值
    n = len(points)
    print("len(points)="+str(n))
    result = points[0][1]
    for i in range(1, n):
        denominator = 1
        for j in range(i):
            denominator= newton_basis(i)
        result += diff_quotient(0,i)*denominator
    return result



# 示例
#points = [(0, 0), (1, 1), (2, 4), (3, 9)]
print("请输入插值节点的个数：")
n = int(input())
points=[]


for i in range(n):
   x = float(input("请输入插值节点的x值："))
   y=float(input("请输入插值节点的y值："))
   
   points.append((x,y))



while True:
    #输出插值结点
    #print("插值节点为：points = [(0, 0), (1, 1), (2, 4), (3, 9)]") 
    print("插值节点为：")
    for i in range(n):
        print("("+str(points[i][0])+","+str(points[i][1])+")")
    print("请输入要计算近似值的点x_interpolate：")
    x_interpolate = float(input())
    result=newton_interpolation(points, x_interpolate)
    print("f("+str(x_interpolate)+")="+str(result))

成效

请输入插值节点的个数：
4
请输入插值节点的x值：1
请输入插值节点的y值：3
请输入插值节点的x值：5
请输入插值节点的y值：5
请输入插值节点的x值：6
请输入插值节点的y值：6
请输入插值节点的x值：7
请输入插值节点的y值：8
插值节点为：
(1.0,3.0)
(5.0,5.0)
(6.0,6.0)
(7.0,8.0)
请输入要计算近似值的点x_interpolate：
3
len(points)=4
f(3.0)=4.4

3.线性函数拟合

def calculate(points):

   def x():
      n=len(points)
      result=0
      for i in range(0,n):
         result +=points[i][0]
      return result
   
   def x2():
      n=len(points)
      result=0
      for i in range(0,n):
        result +=points[i][0]**2
      return result

   def y():
      n=len(points)
      result=0
      for i in range(0,n):
         result +=points[i][1]
      return result
   
   def xy():
      n=len(points)
      result=0
      for i in range(0,n):
         result +=points[i][0]*points[i][1]
      return result
   
   #拟合函数 y=a+b*x
   #a*n+b*x()=y()
   #a*x()+b*x2()=xy()
   #化为
   #a*n*x2()+b*x()*x2()=y()*x2()
   #a*x()*x()+b*x2()*x()=xy()*x()
   n=len(points)
   dividend=n*x2()-x()*x()
   if(dividend==0):
      print("false!float division by zero")
      return 0,0
   
   a=(y()*x2()-x()*xy())/dividend
   b=(y()-a*n)/x()
   return a,b



while True:
   # 示例
   #points = [(0, 0), (1, 1), (2, 4), (3, 9)]
   print("请输入节点的个数：")
   n = int(input())
   points=[]


   for i in range(n):
      x = float(input("请输入节点的x值："))
      y=float(input("请输入节点的y值："))
      
      points.append((x,y))

   #输出插值结点
   #print("插值节点为：points = [(0, 0), (1, 1), (2, 4), (3, 9)]") 
   print("节点为：")
   for i in range(n):
      print("("+str(points[i][0])+","+str(points[i][1])+")")
   #返回拟合函数的a,b的值
   a,b=calculate(points)
   print("拟合函数的a,b的值为："+str(a)+","+str(b))
   print("y="+str(a)+"+"+str(b)+"x")

成效

请输入节点的个数：
4
请输入节点的x值：1
请输入节点的y值：6 
请输入节点的x值：2
请输入节点的y值：7
请输入节点的x值：3
请输入节点的y值：8
请输入节点的x值：4
请输入节点的y值：9
节点为：
(1.0,6.0)
(2.0,7.0)
(3.0,8.0)
(4.0,9.0)
拟合函数的a,b的值为：5.0,1.0
y=5.0+1.0x

收尾

麻了，做完实验感觉还是不太懂，估计只能靠死背公式过及格线了