6-6 Cubic Spline (50分)

最新推荐文章于 2024-07-24 21:32:18 发布

HGGshiwo

最新推荐文章于 2024-07-24 21:32:18 发布

阅读量1.7k

点赞数 3

分类专栏：数值分析

本文链接：https://blog.csdn.net/HGGshiwo/article/details/109678757

版权

数值分析专栏收录该内容

19 篇文章 20 订阅

订阅专栏

Construct the cubic spline interpolant S for the function f, defined at points x0 <x1 <⋯<xn , satisfying some given boundary conditions. Partition a given interval into m equal-length subintervals, and approximate the function values at the endpoints of these subintervals.

Format of functions:

void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, double a[], double b[], double c[], double d[]);
double S( double t, double Fmax, int n, double x[], double a[], double b[], double c[], double d[] );

The function Cubic_Spline is for calculating the set of coefficients of S(x) where S(x)=S[j] (x)=aj +bj (x−xj−1 )+cj (x−xj−1 )2 +dj (x−xj−1 )3 for x∈[xj−1 ,xj ]. The function S is for obtaining an approximation of f(t) using the spline function S(t). The parameters are defined as the following:
int n is the number of interpolating points−1;
double x[] contains n+1 distinct real numbers x0 <x1 <⋯<xn ;
double f[] contains n+1 real numbers f(x0 ),f(x1 ),⋯f(xn );
int Type is the type of boundary conditions, where 1 corresponds to the clamped boundary condition S′(x
0 )=s0 , S′(xn )=sn and 2 corresponds to the natural boundary condition S′′(x0 )=s0 , S′′(xn )=sn ;
double s0 and double sn are s0 and sn , respectively;
double a[], b[], c[], and d[] contain the set of coefficients of S(x);
double t is an endpoint of a subinterval;
and finally double Fmax is the default value that S(t) will assume if t is out of the range [x0 ,xn ].
Sample program of judge:

#include <stdio.h>

#define MAX_N 10

void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, double a[], double b[], double c[], double d[]);

double S( double t, double Fmax, int n, double x[], double a[], double b[], double c[], double d[] );

int main()
{
    int n, Type, m, i;
    double x[MAX_N], f[MAX_N], a[MAX_N], b[MAX_N], c[MAX_N], d[MAX_N];
    double s0, sn, Fmax, t0, tm, h, t;

    scanf("%d", &n);
    for (i=0; i<=n; i++) 
        scanf("%lf", &x[i]);
    for (i=0; i<=n; i++) 
        scanf("%lf", &f[i]);
    scanf("%d %lf %lf %lf", &Type, &s0, &sn, &Fmax);

    Cubic_Spline(n, x, f, Type, s0, sn, a, b, c, d);
    for (i=1; i<=n; i++)
        printf("%12.8e %12.8e %12.8e %12.8e \n", a[i], b[i], c[i], d[i]);

    scanf("%lf %lf %d", &t0, &tm, &m);
    h = (tm-t0)/(double)m;
    for (i=0; i<=m; i++) {
        t = t0+h*(double)i;
        printf("f(%12.8e) = %12.8e\n", t, S(t, Fmax, n, x, a, b, c, d));
    }

    return 0;
}

/* Your functions will be put here */

Sample Input:
2
0.0 1.0 2.0
0.0 1.0 2.0
1 1.0 1.0 0.0
0.0 3.0 2
Sample Output:
0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00
f(0.00000000e+00) = 0.00000000e+00
f(1.50000000e+00) = 1.50000000e+00
f(3.00000000e+00) = 0.00000000e+00

解法：

在这里插入图片描述

1.注意题目中的a,b,c,d下标比书上的都大了1，题目中的范围是1到n，书上的范围是0到n-1，所以下标要调一下。另外a是n+1项，比bcd都多了一项

2.注意书上的算法是s0=0,sn=0所以要修正一下zn,cn

3.如果按照i++去找区间，注意左端点要能算出来。

网上的测试数据：

Sample Input
2

0.0 1.0 2.0

1 1.0 1.0 0.0

0.0 3.0 2

-0.5 -0.25 0.0

-0.0247500 0.3349375 1.1010000

2 0.0 0.0 0.0

-1.0 0.0 4

Sample Output
0.00000000e+000 1.00000000e+000 0.00000000e+000 0.00000000e+000

1.00000000e+000 1.00000000e+000 0.00000000e+000 0.00000000e+000

f(0.00000000e+000) = 0.00000000e+000

f(1.50000000e+000) = 1.50000000e+000

f(3.00000000e+000) = 0.00000000e+000

-2.47500000e-002 1.03237500e+000 0.00000000e+000 6.50200000e+000

3.34937500e-001 2.25150000e+000 4.87650000e+000 -6.50200000e+000

f(-1.00000000e+000) = 0.00000000e+000

f(-7.50000000e-001) = 0.00000000e+000

f(-5.00000000e-001) = -2.47500000e-002

f(-2.50000000e-001) = 3.34937500e-001

f(0.00000000e+000) = 1.10100000e+000

上代码:

#include <stdio.h>

#define MAX_N 10

void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, double a[], double b[], double c[], double d[]);

double S(double t, double Fmax, int n, double x[], double a[], double b[], double c[], double d[]);

int main()
{
    int n, Type, m, i;
    double x[MAX_N], f[MAX_N], a[MAX_N], b[MAX_N], c[MAX_N], d[MAX_N];
    double s0, sn, Fmax, t0, tm, h, t;

    scanf("%d", &n);
    for (i = 0; i <= n; i++)
        scanf("%lf", &x[i]);
    for (i = 0; i <= n; i++)
        scanf("%lf", &f[i]);
    scanf("%d %lf %lf %lf", &Type, &s0, &sn, &Fmax);

    Cubic_Spline(n, x, f, Type, s0, sn, a, b, c, d);
    for (i = 1; i <= n; i++)
        printf("%12.8e %12.8e %12.8e %12.8e \n", a[i], b[i], c[i], d[i]);

    scanf("%lf %lf %d", &t0, &tm, &m);
    h = (tm - t0) / (double)m;
    for (i = 0; i <= m; i++)
    {
        t = t0 + h * (double)i;
        printf("f(%12.8e) = %12.8e\n", t, S(t, Fmax, n, x, a, b, c, d));
    }

    return 0;
}

/* Your functions will be put here */
void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, double a[], double b[], double c[], double d[])
{
    double h[MAX_N];
    double al[MAX_N];
    double l[MAX_N];
    double z[MAX_N];
    double u[MAX_N];

    if (Type == 1)
    { 
        for (int i = 0; i < n; i++)
            h[i] = x[i + 1] - x[i];

        for (int i = 0; i <= n; i++)
            a[i+1] = f[i];
        
        al[0] = 3 * (f[1] - f[0]) / h[0] - 3 * s0;
        al[n] = 3 * sn - 3 * (f[n] - f[n - 1]) / h[n - 1];
        
        for (int i = 1; i < n; i++)
            al[i] = 3 / h[i] * (f[i + 1] - f[i]) - 3 / h[i - 1] * (f[i] - f[i - 1]);

            l[0] = 2 * h[0];
            u[0] = 0.5;
            z[0] = al[0] / l[0];

        for (int i = 1; i < n; i++)
        {
            l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * u[i - 1];
            u[i] = h[i] / l[i];
            z[i] = (al[i] - h[i - 1] * z[i - 1]) / l[i];
        }

        l[n] = h[n - 1] * (2 - u[n - 1]);
        z[n] = (al[n] - h[n - 1] * z[n - 1]) / l[n];
        c[n+1] = z[n];

        for (int j = n - 1; j >= 0; j--)
        {
            c[j+1] = z[j] - u[j] * c[j + 2];
            b[j+1] = (a[j + 2] - a[j+1]) / h[j] - h[j] * (c[j + 2] + 2 * c[j+1]) / 3;
            d[j+1] = (c[j + 2] - c[j+1]) / (3 * h[j]);
        }
    }
    else 
    {   
        for (int i = 0; i < n; i++)
            h[i] = x[i + 1] - x[i];

        for (int i = 0; i <= n; i++)
            a[i+1] = f[i];

        for (int i = 1; i < n; i++)
            al[i] = 3 / h[i] * (f[i + 1] - f[i]) - 3 / h[i - 1] * (f[i] - f[i - 1]);

        l[0] = 1;
        u[0] = 0;
        z[0] = s0/2;
        
        for (int i = 1; i < n; i++)
        {
            l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * u[i - 1];
        u[i] = h[i] / l[i];
        z[i] = (al[i] - h[i - 1] * z[i - 1]) / l[i];
        }

        l[n] = 1;
        z[n] = sn/2;
        c[n+1] = sn/2;

        for (int j = n - 1; j >= 0; j--)
        {
            c[j+1] = z[j] - u[j] * c[j + 2];
            b[j+1] = (a[j + 2] - a[j+1]) / h[j] - h[j]*(c[j + 2] + 2 * c[j+1]) / 3;
            d[j+1] = (c[j + 2] - c[j+1]) / (3 * h[j]);
        }
    }   
}

double S(double t, double Fmax, int n, double x[], double a[], double b[], double c[], double d[])
{
    int i = 0;
    if (x[0]>t||x[n]<t) return Fmax;
    while (t > x[i]) i++;
    if (t == x[0]) i++;
    return a[i] + b[i] * (t - x[i - 1]) + c[i] * (t - x[i - 1]) * (t - x[i - 1]) + d[i] * (t - x[i - 1]) * (t - x[i - 1]) * (t - x[i - 1]);
}