1.txt

关于三次样条插值法求函数的源程序

#include <iostream>
using namespace std;

const int MAXN = 100;

int n;
double x[MAXN], y[MAXN]; //下标从0..n
double alph[MAXN], beta[MAXN], a[MAXN], b[MAXN];
double h[MAXN];
double m[MAXN]; //各点的一阶导数;

inline double sqr(double pa) {
    return pa * pa;
}

double sunc(double p, int i) {
    return (1 + 2 * (p - x[i]) / (x[i + 1] - x[i])) * sqr((p - x[i + 1]) / (x[i + 1] - x[i])) * y[i]
            + (1 + 2 * (p - x[i + 1]) / (x[i] - x[i + 1])) * sqr((p - x[i]) / (x[i + 1] - x[i])) * y[i + 1]
            + (p - x[i]) * sqr((p - x[i + 1]) / (x[i] - x[i + 1])) * m[i]
            + (p - x[i + 1]) * sqr((p - x[i]) / (x[i + 1] - x[i])) * m[i + 1];
}

int main() {
    int i, j;
    double xx;
    freopen("threeInsert.in", "r", stdin);
    scanf("%d", &n);
    for (i = 0; i <= n; i++) scanf("%lf%lf", &x[i], &y[i]);
    // scanf("%lf%lf", &m[0], &m[n]);
    for (i = 0; i <= n - 1; i++) h[i] = x[i + 1] - x[i];
    //第一种边界条件
    //alph[0] = 0; alph[n] = 1; beta[0] = 2 * m[0]; beta[n] = 2 * m[n];
    //第二种边界条件
    alph[0] = 1; alph[n] = 0; beta[0] = 3 * (y[1] - y[0]) / h[0]; beta[n] = 3 * (y[n] - y[n - 1] / h[n - 1]);
    for (i = 1; i <= n - 1; i++) {
        alph[i] = h[i - 1] / (h[i - 1] + h[i]);
        beta[i] = 3 * ((1 - alph[i]) * (y[i] - y[i - 1]) / h[i - 1] + alph[i] * (y[i + 1] - y[i]) / h[i]);
    }
    a[0] = - alph[0] / 2; b[0] = beta[0] / 2;
    for (i = 1; i <= n; i++) {
        a[i] = - alph[i] / (2 + (1 - alph[i]) * a[i - 1]);
        b[i] = (beta[i] - (1 - alph[i]) * b[i - 1]) / (2 + (1 - alph[i]) * a[i - 1]);
    }
    m[n + 1] = 0;
    for (i = n; i >= 0; i--) {
        m[i] = a[i] * m[i + 1] + b[i];
    }
    scanf("%lf", &xx);
    for (i = 0; i < n; i++) {
        if (xx >= x[i] && xx <= x[i + 1]) break;
    }
    printf("%lf\n", sunc(xx, i));
    return 0;
}