📜  贝塞尔插值

📅  最后修改于: 2021-05-06 08:15:17             🧑  作者: Mango

插值是一种针对自变量的任何中间值估算函数值的技术,而计算给定范围之外的函数值的过程称为外推

中心差:中心差运算符d由以下关系定义:

同样,高阶中心差定义为:

注–同一条水平线上的中心差具有相同的后缀

贝塞尔的插值公式–

u = 1/2时,它非常有用。当1/4 时可以给出更好的估计
这里的f(0)是原点,通常取为中点,因为贝塞尔(Bessel’s)用于在中心附近进行插值。
h称为差的间隔,而u =(x – f(0))/ h,其中f(0)是所选原点的项。

例子 –

输入:27.4处的值?

输出 :

27.4的值为3.64968

贝塞尔插值的实现–

C++
// CPP Program to interpolate using Bessel's interpolation
#include 
using namespace std;
  
// calculating u mentioned in the formula
float ucal(float u, int n)
{
    if (n == 0)
        return 1;
  
    float temp = u;
    for (int i = 1; i <= n / 2; i++)
        temp = temp * (u - i);
  
    for (int i = 1; i < n / 2; i++)
        temp = temp * (u + i);
  
    return temp;
}
  
// calculating factorial of given number n
int fact(int n)
{
    int f = 1;
    for (int i = 2; i <= n; i++)
        f *= i;
  
    return f;
}
  
int main()
{
    // Number of values given
    int n = 6;
    float x[] = { 25, 26, 27, 28, 29, 30 };
  
    // y[][] is used for difference table
    // with y[][0] used for input
    float y[n][n];
    y[0][0] = 4.000;
    y[1][0] = 3.846;
    y[2][0] = 3.704;
    y[3][0] = 3.571;
    y[4][0] = 3.448;
    y[5][0] = 3.333;
  
    // Calculating the central difference table
    for (int i = 1; i < n; i++)
        for (int j = 0; j < n - i; j++)
            y[j][i] = y[j + 1][i - 1] - y[j][i - 1];
  
    // Displaying the central difference table
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n - i; j++)
            cout << setw(4) << y[i][j] << "\t";
        cout << endl;
    }
  
    // value to interpolate at
    float value = 27.4;
  
    // Initializing u and sum
    float sum = (y[2][0] + y[3][0]) / 2;
  
    // k is origin thats is f(0)
    int k;
    if (n % 2) // origin for odd
        k = n / 2;
    else
        k = n / 2 - 1; // origin for even
  
    float u = (value - x[k]) / (x[1] - x[0]);
  
    // Solving using bessel's formula
    for (int i = 1; i < n; i++) {
        if (i % 2)
            sum = sum + ((u - 0.5) * 
                  ucal(u, i - 1) * y[k][i]) / fact(i);
        else
            sum = sum + (ucal(u, i) * 
                  (y[k][i] + y[--k][i]) / (fact(i) * 2));
    }
  
    cout << "Value at " << value << " is " << sum << endl;
  
    return 0;
}


Java
// Java Program to interpolate using Bessel's interpolation 
import java.text.*;
class GFG{
// calculating u mentioned in the formula 
static double ucal(double u, int n) 
{ 
    if (n == 0) 
        return 1; 
  
    double temp = u; 
    for (int i = 1; i <= n / 2; i++) 
        temp = temp * (u - i); 
  
    for (int i = 1; i < n / 2; i++) 
        temp = temp * (u + i); 
  
    return temp; 
} 
  
// calculating factorial of given number n 
static int fact(int n) 
{ 
    int f = 1; 
    for (int i = 2; i <= n; i++) 
        f *= i; 
  
    return f; 
} 
  
public static void main(String[] args) 
{ 
    // Number of values given 
    int n = 6; 
    double x[] = { 25, 26, 27, 28, 29, 30 }; 
  
    // y[][] is used for difference table 
    // with y[][0] used for input 
    double[][] y=new double[n][n]; 
    y[0][0] = 4.000; 
    y[1][0] = 3.846; 
    y[2][0] = 3.704; 
    y[3][0] = 3.571; 
    y[4][0] = 3.448; 
    y[5][0] = 3.333; 
  
    // Calculating the central difference table 
    for (int i = 1; i < n; i++) 
        for (int j = 0; j < n - i; j++) 
            y[j][i] = y[j + 1][i - 1] - y[j][i - 1]; 
  
    // Displaying the central difference table
    DecimalFormat df = new DecimalFormat("#.########");
    for (int i = 0; i < n; i++) { 
        for (int j = 0; j < n - i; j++) 
            System.out.print(y[i][j]+"\t"); 
        System.out.println(""); 
    } 
  
    // value to interpolate at 
    double value = 27.4; 
  
    // Initializing u and sum 
    double sum = (y[2][0] + y[3][0]) / 2; 
  
    // k is origin thats is f(0) 
    int k; 
    if ((n % 2)>0) // origin for odd 
        k = n / 2; 
    else
        k = n / 2 - 1; // origin for even 
  
    double u = (value - x[k]) / (x[1] - x[0]); 
  
    // Solving using bessel's formula 
    for (int i = 1; i < n; i++) { 
        if ((i % 2)>0) 
            sum = sum + ((u - 0.5) * 
                ucal(u, i - 1) * y[k][i]) / fact(i); 
        else
            sum = sum + (ucal(u, i) * 
                (y[k][i] + y[--k][i]) / (fact(i) * 2)); 
    } 
  
    System.out.printf("Value at "+value+" is %.5f",sum); 
  
} 
}
// This code is contributed by mits


Python3
# Python3 Program to interpolate
# using Bessel's interpolation
  
# calculating u mentioned in the 
# formula
def ucal(u, n):
  
    if (n == 0):
        return 1;
  
    temp = u;
    for i in range(1, int(n / 2 + 1)):
        temp = temp * (u - i);
  
    for i in range(1, int(n / 2)):
        temp = temp * (u + i);
  
    return temp;
  
# calculating factorial of 
# given number n
def fact(n):
  
    f = 1;
    for i in range(2, n + 1):
        f *= i;
  
    return f;
  
# Number of values given
n = 6;
x = [25, 26, 27, 28, 29, 30];
  
# y[][] is used for difference 
# table with y[][0] used for input
y = [[0 for i in range(n)] 
        for j in range(n)];
y[0][0] = 4.000;
y[1][0] = 3.846;
y[2][0] = 3.704;
y[3][0] = 3.571;
y[4][0] = 3.448;
y[5][0] = 3.333;
  
# Calculating the central
# difference table
for i in range(1, n):
    for j in range(n - i):
        y[j][i] = y[j + 1][i - 1] - y[j][i - 1];
  
# Displaying the central
# difference table
for i in range(n): 
    for j in range(n - i):
        print(y[i][j], "\t", end = " ");
    print("");
  
# value to interpolate at
value = 27.4;
  
# Initializing u and sum
sum = (y[2][0] + y[3][0]) / 2;
  
# k is origin thats is f(0)
k = 0;
if ((n % 2) > 0): # origin for odd
    k = int(n / 2);
else:
    k = int(n / 2 - 1); # origin for even
  
u = (value - x[k]) / (x[1] - x[0]);
  
# Solving using bessel's formula
for i in range(1, n): 
  
    if (i % 2):
        sum = sum + ((u - 0.5) *
                 ucal(u, i - 1) *
              y[k][i]) / fact(i);
    else:
        sum = sum + (ucal(u, i) * (y[k][i] + 
                     y[k - 1][i]) / (fact(i) * 2));
        k -= 1;
  
print("Value at", value, "is", round(sum, 5));
  
# This code is contributed by mits


C#
// C# Program to interpolate using Bessel's interpolation 
  
class GFG{
// calculating u mentioned in the formula 
static double ucal(double u, int n) 
{ 
    if (n == 0) 
        return 1; 
  
    double temp = u; 
    for (int i = 1; i <= n / 2; i++) 
        temp = temp * (u - i); 
  
    for (int i = 1; i < n / 2; i++) 
        temp = temp * (u + i); 
  
    return temp; 
} 
  
// calculating factorial of given number n 
static int fact(int n) 
{ 
    int f = 1; 
    for (int i = 2; i <= n; i++) 
        f *= i; 
  
    return f; 
} 
  
public static void Main() 
{ 
    // Number of values given 
    int n = 6; 
    double []x = { 25, 26, 27, 28, 29, 30 }; 
  
    // y[,] is used for difference table 
    // with y[,0] used for input 
    double[,] y=new double[n,n]; 
    y[0,0] = 4.000; 
    y[1,0] = 3.846; 
    y[2,0] = 3.704; 
    y[3,0] = 3.571; 
    y[4,0] = 3.448; 
    y[5,0] = 3.333; 
  
    // Calculating the central difference table 
    for (int i = 1; i < n; i++) 
        for (int j = 0; j < n - i; j++) 
            y[j,i] = y[j + 1,i - 1] - y[j,i - 1]; 
  
    // Displaying the central difference table
    for (int i = 0; i < n; i++) { 
        for (int j = 0; j < n - i; j++) 
            System.Console.Write(y[i,j]+"\t"); 
        System.Console.WriteLine(""); 
    } 
  
    // value to interpolate at 
    double value = 27.4; 
  
    // Initializing u and sum 
    double sum = (y[2,0] + y[3,0]) / 2; 
  
    // k is origin thats is f(0) 
    int k; 
    if ((n % 2)>0) // origin for odd 
        k = n / 2; 
    else
        k = n / 2 - 1; // origin for even 
  
    double u = (value - x[k]) / (x[1] - x[0]); 
  
    // Solving using bessel's formula 
    for (int i = 1; i < n; i++) { 
        if ((i % 2)>0) 
            sum = sum + ((u - 0.5) * 
                ucal(u, i - 1) * y[k,i]) / fact(i); 
        else
            sum = sum + (ucal(u, i) * 
                (y[k,i] + y[--k,i]) / (fact(i) * 2)); 
    } 
  
    System.Console.WriteLine("Value at "+value+" is "+System.Math.Round(sum,5)); 
  
} 
}
// This code is contributed by mits


PHP


输出:

4    -0.154    0.0120001    -0.00300002    0.00399971    -0.00699902    
3.846    -0.142    0.00900006    0.000999689    -0.00299931    
3.704    -0.133    0.00999975    -0.00199962    
3.571    -0.123    0.00800014    
3.448    -0.115    
3.333    
Value at 27.4 is 3.64968