📜  斯特林插值公式程序

📅  最后修改于: 2021-04-30 02:37:13             🧑  作者: Mango

给定n个浮点值x及其对应的函数值f(x),针对自变量x的任意中间值(即,在x = a时)估算数学函数的值。
例子:

输入:n = 5 x _1 = 0,x _2 = 0.5,x _3 = 1.0,x _4 = 1.5,x _5 = 2.0 f(x _1 )= 0,f(x _2 )= 0.191,f(x _3 )= 0.341,f(x _4 )= 0.433,f(x _5 )= 0.477 a = 1.22输出:函数1.22的值为0.389。可以看出f(1.0)= 0.341和f(1.5)= 0.433,因此f(1.22)应该介于这两个值之间。使用斯特林近似,f(1.22)为0.389。输入:n = 7 x _1 = 0,x _2 = 5,x _3 = 10,x _4 = 15,x _5 = 20,x _6 = 25,x _7 = 30 f(x _1 )= 0,f(x _2 )= 0.0875,f(x _3 )= 0.1763,f(x _4 )= 0.2679,f(x _5 )= 0.364,f(x _6 )= 0.4663,f(x _7 )= 0.5774 a = 16输出:16处的函数值为0.2866。

斯特林实习

斯特林近似或斯特林插值公式是一种插值技术,用于获取已知数据点的离散集合范围内的中间点处的函数值。

斯特林公式是通过采用高斯正向和高斯反向公式的平均值或均值来获得的。高斯前进和后退公式都是用于在列表集的中部附近获得函数值的公式。

怎么找

斯特林近似法涉及使用前向差分表,该表可以从给定的x和f(x)或y的给定集合中准备,如下所示–

该表是在x及其对应的f(x)或y的帮助下准备的。然后,通过计算上一列中其前一个值和后一个值之间的差来计算每个下一个列的值,例如\Delta y_0 = y _1 – y _0\Delta y_1 = y _2 – y _1\Delta ^2 y_0 = \Delta y_1\Delta y_0 , 等等。

现在,用于在a处获得f(x)或y的高斯正向公式为:
  y_p = y_0 + p\Delta y_0 + \frac{p(p-1)}{2!}\Delta ^2 y_-_1+\frac{p(p-1)(p+1)}{3!}\Delta ^3 y_-_1 + \frac{p(p-1)(p+1)(p-2)}{4!}\Delta ^4 y_-_2 + .......

在哪里,
p = \frac{a - x_0}{h}
a是我们必须确定f(x),x的点_0是从给定的x中更接近a的选定值(通常是从表的中间选择一个值),h是任意两个连续x之差。现在,y _0变成与x对应的值_0和x之前的值_0下标为负,之后的为x_0下标为正,如下表所示–

用于在a处获得f(x)或y的高斯向后公式为:

  y_p = y_0 + p\Delta y_-_1 + \frac{p(p+1)}{2!}\Delta ^2 y_-_1+\frac{p(p-1)(p+1)}{3!}\Delta ^3 y_-_2 + \frac{p(p+1)(p-1)(p+2)}{4!}\Delta ^4 y_-_2 + .......

现在,以上述两个公式的平均值为基础,并获得斯特林近似公式,如下所示:

  P_2_n(x) = y_0 + q.\frac{\Delta y_-_1+\Delta y_0}{2}+\frac{q^2}{2!},\Delta^2 y_-_1 +\frac{q(q^2-1)}{3!}.\frac{\Delta^3 y_-_2 + \Delta^3 y_-_1}{2}+\frac{q^2(q^2-1)}{4!}. \Delta ^4 y_-_2 + \frac{q(q^2-1)(q^2-2^2)}{5!}.\frac{\Delta ^5 y_-_3+ \Delta ^5y_-_2}{2}+\frac{q^2(q^2-1)(q^2-2^2)}{6!}.\Delta^6y_-_3+....+ \frac{q(q^2-1)(q^2-2^2)(q^2-3^2)...[q^2-(n-1)^2]}{(2n-1)!}* \frac{\Delta ^{2n-1}y_-_n+ \Delta ^{2n-1}y_-_{n-1}}{2}+\frac{q^2(q^2-1)(q^2-2^2)...[q^2-(n-1)^2]}{(2n)!} \Delta^{2n}y_{-n},

此处,q与高斯公式中的p相同,其余所有符号均相同。

何时使用

  • 当q介于\frac{-1}{2}\frac{1}{2}.
  • 超出此范围,仍可以使用,但计算值的准确性会降低。
  • 它给出了范围的最佳估计\frac{-1}{4}
C++
// C++ program to demonstrate Stirling
// Approximation
#include 
using namespace std;
  
// Function to calculate value using 
// Stirling formula
void Stirling(float x[], float fx[], float x1,
                                    int n)
{
    float h, a, u, y1 = 0, N1 = 1, d = 1,
    N2 = 1, d2 = 1, temp1 = 1, temp2 = 1,
    k = 1, l = 1, delta[n][n];
    
    int i, j, s;
    h = x[1] - x[0];
    s = floor(n / 2);
    a = x[s];
    u = (x1 - a) / h;
  
    // Preparing the forward difference
    // table
    for (i = 0; i < n - 1; ++i) {
        delta[i][0] = fx[i + 1] - fx[i];
    }
    for (i = 1; i < n - 1; ++i) {
        for (j = 0; j < n - i - 1; ++j) {
            delta[j][i] = delta[j + 1][i - 1]
                          - delta[j][i - 1];
        }
    }
  
    // Calculating f(x) using the Stirling 
    // formula
    y1 = fx[s];
  
    for (i = 1; i <= n - 1; ++i) {
        if (i % 2 != 0) {
            if (k != 2) {
                temp1 *= (pow(u, k) - 
                          pow((k - 1), 2));
            } 
            else {
                temp1 *= (pow(u, 2) - 
                          pow((k - 1), 2));
            }
            ++k;
            d *= i;
            s = floor((n - i) / 2);
            y1 += (temp1 / (2 * d)) * 
                   (delta[s][i - 1] + 
                   delta[s - 1][i - 1]);
        } 
        else {
            temp2 *= (pow(u, 2) - 
                      pow((l - 1), 2));
            ++l;
            d *= i;
            s = floor((n - i) / 2);
            y1 += (temp2 / (d)) *
                  (delta[s][i - 1]);
        }
    }
  
    cout << y1;
}
  
// Driver main function
int main()
{
    int n;
    n = 5;
    float x[] = { 0, 0.5, 1.0, 1.5, 2.0 };
    float fx[] = { 0, 0.191, 0.341, 0.433,
                             0.477 };
  
    // Value to calculate f(x)
    float x1 = 1.22;
  
    Stirling(x, fx, x1, n);
    return 0;
}


Java
// Java program to demonstrate Stirling 
// Approximation
import java.io.*;
import static java.lang.Math.*;
  
public class A {
    static void Stirling(double x[], double fx[],
                         double x1, int n)
    {
        double h, a, u, y1 = 0, N1 = 1, d = 1,
            N2 = 1, d2 = 1, temp1 = 1,
          temp2 = 1, k = 1, l = 1, delta[][];
          
        delta = new double[n][n];
        int i, j, s;
        h = x[1] - x[0];
        s = (int)floor(n / 2);
        a = x[s];
        u = (x1 - a) / h;
  
        // Preparing the forward difference 
       // table
        for (i = 0; i < n - 1; ++i) {
            delta[i][0] = fx[i + 1] - fx[i];
        }
        for (i = 1; i < n - 1; ++i) {
            for (j = 0; j < n - i - 1; ++j) {
                delta[j][i] = delta[j + 1][i - 1] 
                              - delta[j][i - 1];
            }
        }
  
        // Calculating f(x) using the Stirling 
        // formula
        y1 = fx[s];
  
        for (i = 1; i <= n - 1; ++i) {
            if (i % 2 != 0) {
                if (k != 2) {
                    temp1 *= (pow(u, k) - 
                              pow((k - 1), 2));
                } 
                else {
                    temp1 *= (pow(u, 2) - 
                              pow((k - 1), 2));
                }
                ++k;
                d *= i;
                s = (int)floor((n - i) / 2);
                y1 += (temp1 / (2 * d)) * 
                     (delta[s][i - 1] +
                      delta[s - 1][i - 1]);
            } 
            else {
                temp2 *= (pow(u, 2) -
                        pow((l - 1), 2));
                ++l;
                d *= i;
                s = (int)floor((n - i) / 2);
                y1 += (temp2 / (d)) *
                      (delta[s][i - 1]);
            }
        }
  
        System.out.print(+ y1);
    }
  
    // Driver main function
public static void main(String args[])
    {
        int n;
        n = 5;
        double x[] = {0, 0.5, 1.0, 1.5, 2.0 };
        double fx[] = {0, 0.191, 0.341, 0.433,
                                      0.477 };
  
        // Value to calculate f(x)
        double x1 = 1.22;
  
        Stirling(x, fx, x1, n);
    }
}


Python3
# Python3 program to demonstrate Stirling
# Approximation
import math
  
# Function to calculate value using 
# Stirling formula
def Stirling(x, fx, x1, n):
  
    y1 = 0; N1 = 1; d = 1;
    N2 = 1; d2 = 1; 
    temp1 = 1; temp2 = 1;
    k = 1; l = 1; 
    delta = [[0 for i in range(n)]
                for j in range(n)];
  
    h = x[1] - x[0];
    s = math.floor(n / 2);
    a = x[s];
    u = (x1 - a) / h;
  
    # Preparing the forward difference
    # table
    for i in range(n - 1): 
        delta[i][0] = fx[i + 1] - fx[i];
    for i in range(1, n - 1):
        for j in range(n - i - 1):
            delta[j][i] = (delta[j + 1][i - 1] - 
                           delta[j][i - 1]);
  
    # Calculating f(x) using the Stirling formula
    y1 = fx[s];
  
    for i in range(1, n):
        if (i % 2 != 0): 
            if (k != 2): 
                temp1 *= (pow(u, k) - pow((k - 1), 2));
            else:
                temp1 *= (pow(u, 2) - pow((k - 1), 2));
            k += 1;
            d *= i;
            s = math.floor((n - i) / 2);
            y1 += (temp1 / (2 * d)) * (delta[s][i - 1] + 
                                       delta[s - 1][i - 1]);
        else:
            temp2 *= (pow(u, 2) - pow((l - 1), 2));
            l += 1;
            d *= i;
            s = math.floor((n - i) / 2);
            y1 += (temp2 / (d)) * (delta[s][i - 1]);
  
    print(round(y1, 3));
  
# Driver Code
n = 5;
x = [0, 0.5, 1.0, 1.5, 2.0 ];
fx = [ 0, 0.191, 0.341, 0.433, 0.477];
  
# Value to calculate f(x)
x1 = 1.22;
Stirling(x, fx, x1, n);
  
# This code is contributed by mits


C#
// C# program to demonstrate Stirling 
// Approximation
using System;
  
public class A 
{
static void Stirling(double[] x, double[] fx,
                     double x1, int n)
{
    double h, a, u, y1 = 0, d = 1, temp1 = 1,
                     temp2 = 1, k = 1, l = 1;
    double[,] delta;
      
    delta = new double[n, n];
    int i, j, s;
    h = x[1] - x[0];
    s = (int)Math.Floor((double)(n / 2));
    a = x[s];
    u = (x1 - a) / h;
  
    // Preparing the forward difference 
    // table
    for (i = 0; i < n - 1; ++i) 
    {
        delta[i, 0] = fx[i + 1] - fx[i];
    }
    for (i = 1; i < n - 1; ++i) 
    {
        for (j = 0; j < n - i - 1; ++j) 
        {
            delta[j, i] = delta[j + 1, i - 1] 
                        - delta[j, i - 1];
        }
    }
  
    // Calculating f(x) using the Stirling 
    // formula
    y1 = fx[s];
  
    for (i = 1; i <= n - 1; ++i) 
    {
        if (i % 2 != 0) 
        {
            if (k != 2) 
            {
                temp1 *= (Math.Pow(u, k) - 
                          Math.Pow((k - 1), 2));
            } 
            else 
            {
                temp1 *= (Math.Pow(u, 2) - 
                          Math.Pow((k - 1), 2));
            }
            ++k;
            d *= i;
            s = (int)Math.Floor((double)((n - i) / 2));
            y1 += (temp1 / (2 * d)) * 
                  (delta[s, i - 1] +
                   delta[s - 1, i - 1]);
        } 
        else 
        {
            temp2 *= (Math.Pow(u, 2) -
                      Math.Pow((l - 1), 2));
            ++l;
            d *= i;
            s = (int)Math.Floor((double)((n - i) / 2));
            y1 += (temp2 / (d)) *
                  (delta[s, i - 1]);
        }
    }
  
    Console.Write(+ y1);
}
  
// Driver Code
public static void Main()
{
    int n;
    n = 5;
    double[] x = {0, 0.5, 1.0, 1.5, 2.0 };
    double[] fx = {0, 0.191, 0.341, 0.433,
                                0.477 };
  
    // Value to calculate f(x)
    double x1 = 1.22;
  
    Stirling(x, fx, x1, n);
}
}
  
// This code is contributed 
// by Akanksha Rai


PHP


输出:

0.389 

斯特林公式相对于其他相似公式的主要优点是,与其他差分公式相比,其下降速度要快得多,因此考虑到前几个项本身会提供更高的精度,而它的缺点是,对于斯特林逼近而言,应该有任何两个连续x之间的均等差。

参考-BS Grewal的《高级工程数学》。