三次样条函数三弯矩算法
- 格式：docx
- 大小：56.08 KB
- 文档页数：7
摘要
求函数在给定区间上的定积分，在微积分学中已给出了许多计算方法，但是，在实际问题计算中，往往仅给出函数在一些离散点的值，它的解析表达式没有明显的给出，或者，虽然给出解析表达式，但却很难求得其原函数。这时我们可以通过数值方法求出函数积分的近似值。
在用近似值代替真实值时，遇到的问题就是近似值的代数精度是否足够。当代数精度不足够时，很显然提高插值函数的次数是一种方法，但是考虑到数值计算的稳定性，当次数过高时，会出现龙格现象，用增大n的方法来提高数值积代数精度是不可取的。
正如我们所知道的分段线性插值，逼近程度好，但光滑性差。分段三次Hermite插值，逼近程度好，光滑性也有所提高，但是需要增加更多的条件，不太实用。因此，我们将介绍一种结合二者的优点的插值方法——三次样条插值。本实验将介绍三次样条插值的三弯矩算法。
关键词：龙格现象 三弯矩算法 代数精度 分段三次Hermite插值
1 1、实验目的
1) 通过本次实验体会并学习三次样条插值的优点。
2) 通过对三次样条插值进行编程实现，提高自己的编程能力。
3) 用实验报告的形式展现，提高自己在写论文方面的能力。
2、算法流程
如果已知函数)(xfy在节点a=𝑥0<𝑥1<⋯<𝑥𝑛=𝑏,𝑦𝑖=𝑓(𝑥𝑖),𝑖=0,1,2,⋯,𝑛处的函数值和导数值：
nixfyii,,2,1,0),(
如果)(xS满足条件：
1) S(x)是一个分段的三次多项式且iiyxs)(；
2) S(x)在[a,b]具有二阶连续导数。
则称S(x)是三次样条插值函数。S(x)的具体形式为：
其中S𝑖(𝑥)在[𝑥𝑛−1,𝑥𝑛]上是三次多项式S𝑖(x)=𝑎𝑖𝑥3+𝑏𝑖𝑥2+𝑐𝑖𝑥+𝑑𝑖
由插值条件S(𝑥𝑖)=𝑦𝑖，i=0,1,2,…,n，得n+1个条件。
边界条件一：𝑆′(𝑥0)=𝑦0′,𝑆′(𝑥𝑛)=𝑦𝑛′
边界条件二：𝑆′′(𝑥0)=𝑦0′′,𝑆′′(𝑥𝑛)=𝑦𝑛′′
边界条件三：假定函数𝑦=𝑓(𝑥)是以b-a为周期的周期函数，这时要求S(x)也是周期函数，即
{𝑆(𝑥0+0)=𝑆(𝑥𝑛−0)𝑆′(𝑥0+0)=𝑆′(𝑥𝑛−0)𝑆′′(𝑥0+0)=𝑆′′(𝑥𝑛+0)
]12,121,01,[,...............][,][,nnnxxxxsxxxxsxxxxsxs2 针对三种边界条件的求解方法的不同，可以分为三转角算法和三弯矩算法，本实验将介绍和学习三转角算法。
三弯矩算法：
在每个子区间[𝑥𝑛−1,𝑥𝑛]，因为S(x)是三次多项式，所以S′′(x)是一次多项式，假设节点𝑥𝑖处
S′′(𝑥𝑖)=𝑀𝑖,𝑖=0,1,2,⋯,𝑛
则在[𝑥𝑛−1,𝑥𝑛]上
𝑆𝑖′′(x)=𝑥−𝑥𝑖𝑥𝑖−1−𝑥𝑖𝑀𝑖−1+𝑥−𝑥𝑖−1𝑥𝑖−𝑥𝑖−1𝑀𝑖=𝑥𝑖−𝑥ℎ𝑖𝑀𝑖−1+𝑥−𝑥𝑖−1ℎ𝑖𝑀𝑖
对上式积分两次得
𝑆𝑖(𝑥)=𝑀𝑖−1(𝑥𝑖−𝑥)36ℎ𝑖+𝑀𝑖(𝑥−𝑥𝑖−1)36ℎ𝑖+𝐶𝑖(𝑥𝑖−𝑥)+𝐷𝑖(𝑥−𝑥𝑖−1)
其中𝐶𝑖和𝐷𝑖的值为任意常数，由端点值可以求得。
显然只要确定出𝑀𝑖(𝑖=0,1,2,⋯,𝑛)就能求出𝑆𝑖(𝑥)，进而得到三次样条函数S(x)。关于𝑀𝑖，我们利用样条函数在节点𝑥𝑖处一阶导数连续来确定。
3、算法实例
已知函数𝑦=𝑓(𝑥)的如下表的数据
i 0 1 2 3
𝑥𝑖 1 2 4 5
𝑦𝑖 1 3 4 2
求满足自然边界条件S′′(𝑥0)=S′′(𝑥3)=0的𝑓(3)的近似值
解：根据题意可知，函数在区间两端的二阶导数，因此可知应该用三弯矩算法来实现，具体程序如下：
#include "stdio.h"
#define N 4
void main()
{
int i,k;
float X,s,y0,yn;
float a[N][N+1],h[N],u[N],v[N],g[N],m[N],p[N],q[N],w[N];
printf("please input X:"); //X为未知数的大小
scanf("%f",&X); 3 printf("please input x:"); //输入x的大小
for(i=0;i
scanf("%f",&a[i][0]);
printf("please input y:"); //输入y的大小
for(i=0;i
scanf("%f",&a[i][1]);
for(i=1;i
h[i]=a[i][0]-a[i-1][0]; //计算步长
for(i=1;i
{
v[i]=h[i+1]/(h[i]+h[i+1]);
u[i]=1-v[i];
g[i]=3*u[i]*(a[i+1][1]-a[i][1])/h[i+1]+3*v[i]*(a[i][1]-a[i-1][1])/h[i];
}
printf("\t(1)已知边界条件1\n");
printf("\t(2)已知边界条件2\n");
printf("请选择边界条件序号:");
scanf("%d",&k);
if(k==1)
{
printf("请输入y0和yn的一阶导:"); //输入边界条件一
scanf("%f%f",&m[0],&m[N-1]);
p[0]=0; //用追赶法求解m[N]
q[0]=0;
g[1]=g[1]-v[1]*m[0];
g[N-2]=g[N-2]-u[N-2]*m[N-1];
for(i=1;i
{
w[i]=2-u[i]*p[i-1];
p[i]=v[i]/w[i];
q[i]=(g[i]-u[i]*q[i-1])/w[i];
}
m[N-2]=q[N-2];
for(i=N-3;i>0;i--) 4 m[i]=q[i]-p[i]*m[i+1];
printf("输出m[i]的值:\n");
for(i=0;i
printf("%f\n",m[i]);
for(i=1;i
if(X>a[i-1][0]&&X
s=(h[i]+2*(X-a[i-1][0]))*(X-a[i][0])*(X-a[i][0])*a[i-1][1]/(h[i]*h[i]*h[i])+(h[i]-2*(X-a[i][0]))*(X-a[i-1][0])*(X-a[i-1][0])*a[i][1]/(h[i]*h[i]*h[i])+(X-a[i-1][0])*(X-a[i][0])*(X-a[i][0])*m[i-1]/(h[i]*h[i])+(X-a[i-1][0])*(X-a[i-1][0])*(X-a[i][0])*m[i]/(h[i]*h[i]);
printf("s(%f)=%f\n",X,s);
}
if(k==2)
{
printf("请输入y0和yn的二阶导:"); //输入边界条件二
scanf("%f%f",&y0,&yn);
g[0]=3*(a[1][1]-a[0][1])/h[1]-h[1]*y0/2;
g[N-1]=3*(a[N-1][1]-a[N-2][1])/h[N-1]+h[N-1]*yn/2;
q[0]=g[0];
u[0]=1;
v[N-1]=1;
w[0]=2;
for(i=1;i
{
w[i]=2-v[i]*u[i-1]/w[i-1];
q[i]=g[i]-v[i]*q[i-1]/w[i-1];
}
m[N-1]=q[N-1]/w[N-1];
for(i=N-2;i>=0;i--)
m[i]=(q[i]-u[i]*m[i+1])/w[i];
printf("输出m[i]的值:\n");
for(i=0;i
printf("%f\n",m[i]);
for(i=1;i
if(X>=a[i-1][0]&&X<=a[i][0]) 5
s=(h[i]+2*(X-a[i-1][0]))*(X-a[i][0])*(X-a[i][0])*a[i-1][1]/(h[i]*h[i]*h[i])+(h[i]-2*(X-a[i][0]))*(X-a[i-1][0])*(X-a[i-1][0])*a[i][1]/(h[i]*h[i]*h[i])+(X-a[i-1][0])*(X-a[i][0])*(X-a[i][0])*m[i-1]/(h[i]*h[i])+(X-a[i-1][0])*(X-a[i-1][0])*(X-a[i][0])*m[i]/(h[i]*h[i]);
printf("s(%f)=%f\n",X,s);
}
}
运行结果：
5、对结果进行分析
通过用编程实现对上例的求解，可以看出结果较为准确，而且不仅逼近效果好，而且光滑性有所提高，但是由于在计算机上计算，会存在计算误差。
6、参考文献
[1] 秦新强.数值逼近.西安：西安理工大学出版社，2010
