首页 新闻 会员 周边 捐助

Cubic Spline Interpolation (三次样条插值)

0
悬赏园豆:20 [待解决问题]

代码转载来自: https://gist.github.com/dreikanter/3526685
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;

namespace Interpolation
{
/// <summary>
/// Spline interpolation class.
/// </summary>
public class SplineInterpolator
{
private readonly double[] _keys;

    private readonly double[] _values;
    
    private readonly double[] _h;
    
    private readonly double[] _a;

    /// <summary>
    /// Class constructor.
    /// </summary>
    /// <param name="nodes">Collection of known points for further interpolation.
    /// Should contain at least two items.</param>
    public SplineInterpolator(IDictionary<double, double> nodes)
    {
        if (nodes == null)
        {
            throw new ArgumentNullException("nodes");
        }

        var n = nodes.Count;

        if (n < 2)
        {
            throw new ArgumentException("At least two point required for interpolation.");
        }

        _keys = nodes.Keys.ToArray();
        _values = nodes.Values.ToArray();
        _a = new double[n];
        _h = new double[n];

        for (int i = 1; i < n; i++)
        {
            _h[i] = _keys[i] - _keys[i - 1];
        }

        if (n > 2)
        {
            var sub = new double[n - 1];
            var diag = new double[n - 1];
            var sup = new double[n - 1];

            for (int i = 1; i <= n - 2; i++)
            {
                diag[i] = (_h[i] + _h[i + 1]) / 3;
                sup[i] = _h[i + 1] / 6;
                sub[i] = _h[i] / 6;
                _a[i] = (_values[i + 1] - _values[i]) / _h[i + 1] - (_values[i] - _values[i - 1]) / _h[i];
            }

            SolveTridiag(sub, diag, sup, ref _a, n - 2);
        }
    }

    /// <summary>
    /// Gets interpolated value for specified argument.
    /// </summary>
    /// <param name="key">Argument value for interpolation. Must be within 
    /// the interval bounded by lowest ang highest <see cref="_keys"/> values.</param>
    public double GetValue(double key)
    {
        int gap = 0;
        var previous = double.MinValue;

        // At the end of this iteration, "gap" will contain the index of the interval
        // between two known values, which contains the unknown z, and "previous" will
        // contain the biggest z value among the known samples, left of the unknown z
        for (int i = 0; i < _keys.Length; i++)
        {
            if (_keys[i] < key && _keys[i] > previous)
            {
                previous = _keys[i];
                gap = i + 1;
            }
        }

        var x1 = key - previous;
        var x2 = _h[gap] - x1;

        return ((-_a[gap - 1] / 6 * (x2 + _h[gap]) * x1 + _values[gap - 1]) * x2 +
            (-_a[gap] / 6 * (x1 + _h[gap]) * x2 + _values[gap]) * x1) / _h[gap];
    }


    /// <summary>
    /// Solve linear system with tridiagonal n*n matrix "a"
    /// using Gaussian elimination without pivoting.
    /// </summary>
    private static void SolveTridiag(double[] sub, double[] diag, double[] sup, ref double[] b, int n)
    {
        int i;

        for (i = 2; i <= n; i++)
        {
            sub[i] = sub[i] / diag[i - 1];
            diag[i] = diag[i] - sub[i] * sup[i - 1];
            b[i] = b[i] - sub[i] * b[i - 1];
        }

        b[n] = b[n] / diag[n];
        
        for (i = n - 1; i >= 1; i--)
        {
            b[i] = (b[i] - sup[i] * b[i + 1]) / diag[i];
        }
    }
}

}

想问一下这个代码的这个部分
for (int i = 1; i <= n - 2; i++)
{
diag[i] = (_h[i] + _h[i + 1]) / 3;
sup[i] = _h[i + 1] / 6;
sub[i] = _h[i] / 6;
_a[i] = (_values[i + 1] - _values[i]) / _h[i + 1] - (_values[i] - _values[i - 1]) / _h[i];
}

这个部分我看不懂,想请问一下为什么diagonal需要除3,superdiagonal需要除6,subdiagonal为什么需要除6, 还有_a 这个变量 是拿来做什么的,?
有没有大牛可以帮忙解释一下,还有这个方程式来自哪里?

SolveTridiag()这个方法的代码我看的懂,网上有方程式提供.

接下来是 public double GetValue(double key) 方法的代码, 我也看不懂.

return ((-_a[gap - 1] / 6 * (x2 + _h[gap]) * x1 + _values[gap - 1]) * x2 +
(-_a[gap] / 6 * (x1 + _h[gap]) * x2 + _values[gap]) * x1) / _h[gap];

特别是这一句, 有没有大牛可以解释一下, 顺便告知代码的方程式来自哪里?

我查找了好的关于三次样条插值的资料跟方程式, 可是没有查找到类似的方程式运用在这个代码上.

麻烦大牛帮帮忙解我心中的疑惑!!!

辰洁的主页 辰洁 | 初学一级 | 园豆:182
提问于:2021-08-18 20:14
< >
分享
清除回答草稿
   您需要登录以后才能回答,未注册用户请先注册