📜  高斯前向插值

📅  最后修改于: 2021-09-23 04:39:52             🧑  作者: Mango

插值是指在给定的数据集中创建新数据点的过程。下面的代码使用高斯给出的公式计算给定范围内的离散数据集所需的数据点,这种方法称为高斯前向法。

高斯前向法:

高斯插值属于中心差分插值公式。假设对于一组 x 值,我们给出了以下 y=f(x) 值:
X: x0 x1 x2 ………. xn
Y: y0 y1 y2 ………… yn
分别用Δy0、Δy1、Δy2、……、Δyn-1表示的差y1-y0、y2-y1、y3-y2、……、yn-yn-1,称为一阶前向差。因此,第一个前向差异是:
Δy 0 = y 1 – y 0
以同样的方式我们可以计算高阶差异。

在创建表之后,我们根据以下公式计算值:

现在,让我们举一个例子来解决它,以便更好地理解。
问题:
从下表中,使用高斯正向公式找到e 1.17的值。

x 1.00 1.05 1.10 1.15 1.20 1.25 1.30
ex 2.7183 2.8577 3.0042 3.1582 3.3201 3.4903 3.6693

解决方案:
我们有

y p = y 0 + pΔy 0 + (p(p-1)/2!).Δy 2 0 + ((p+1)p(p-1)/3!).Δy 3 0 + …
其中 p = (x 1.17 – x 1.15 ) / h
并且 h = x 1 – x 0 = 0.05
所以,p = 0.04
现在,我们需要计算 Δy 0 、Δy 2 0 、Δy 3 0 ……等等。

将所需的值放入公式中-
y x = 1.17 = 3.158 + (2/5)(0.162) + (2/5)(2/5 – 1)/2.(0.008) …
y x = 1.17 = 3.2246
代码:实现高斯前向公式的Python代码

Python3
# Python3 code for Gauss's Forward Formula
# importing library
import numpy as np
 
# function for calculating coefficient of Y
def p_cal(p, n):
 
    temp = p;
    for i in range(1, n):
         if(i%2==1):
             temp * (p - i)
         else:
             temp * (p + i)
    return temp;
# function for factorial
def fact(n):
    f = 1
    for i in range(2, n + 1):
        f *= i
    return f
 
# storing available data
n = 7;
x = [ 1, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30 ];
 
y = [[0 for i in range(n)]
        for j in range(n)];
y[0][0] = 2.7183;
y[1][0] = 2.8577;
y[2][0] = 3.0042;
y[3][0] = 3.1582; 
y[4][0] = 3.3201;
y[5][0] = 3.4903;
y[6][0] = 3.6693;
 
# Generating Gauss's triangle
for i in range(1, n):
    for j in range(n - i):
        y[j][i] = np.round((y[j + 1][i - 1] - y[j][i - 1]),4);
 
# Printing the Triangle
for i in range(n):
    print(x[i], end = "\t");
    for j in range(n - i):
        print(y[i][j], end = "\t");
    print("");
 
# Value of Y need to predict on
value = 1.17;
 
# implementing Formula
sum = y[int(n/2)][0];
p = (value - x[int(n/2)]) / (x[1] - x[0])
 
for i in range(1,n):
    # print(y[int((n-i)/2)][i])
    sum = sum + (p_cal(p, i) * y[int((n-i)/2)][i]) / fact(i)
 
print("\nValue at", value,
    "is", round(sum, 4));


输出 :

1       2.7183  0.1394  0.0071  0.0004  0.0     0.0     0.0001  
1.05    2.8577  0.1465  0.0075  0.0004  0.0     0.0001  
1.1     3.0042  0.154   0.0079  0.0004  0.0001  
1.15    3.1582  0.1619  0.0083  0.0005  
1.2     3.3201  0.1702  0.0088  
1.25    3.4903  0.179   
1.3     3.6693  

Value at 1.17 is 3.2246