三次样条函数三弯矩算法

1、摘要求函数在给定区间上的定积分，在微积分学中已给出了许多计算方法, 但是，在实际问题计算中，往往仅给出函数在一些离散点的值，它的解析 表达式没有明显的给出，或者，虽然给出解析表达式，但却很难求得其原 函数，这时我们可以通过数值方法求出函数积分的近似值。在用近似值代替真实值时，遇到的问题就是近似值的代数精度是否足 够。当代数精度不足够时，很显然提高插值函数的次数是一种方法，但是 考虑到数值计算的稳定性，当次数过高时，会出现龙格现象，用增大n的 方法来提高数值积代数精度是不可取的。正如我们所知道的分段线性插值，逼近程度好，但光滑性差。分段三 次Hermite插值，逼近程度好，光滑性也有所提高，但是

2、需要增加更多的 条件，不太实用。因此，我们将介绍一种结合二者的优点的插值方法一一 三次样条插值。本实验将介绍三次样条插值的三弯矩算法。关键词：龙格现象三弯矩算法代数精度分段三次Hermite插值1、实验目的1)通过本次实验体会并学习三次样条插值的优点。2)通过对三次样条插值进行编程实现，提高自己的编程能力。3)用实验报告的形式展现，提高自己在写论文方面的能力。2、算法流程如果已知函数y= f(x)在节点a = %o <<V = b,yt = 0,1,2,，几处的函数值和导数值：y = f("),i = 0,LAm如果S(x)满足条件:1) S(x)是一个分段的三次多项式且

3、*(x) = v ；2) S(x)在a,b具有二阶连续导数。则称S(x)是三次样条插值函数。S(x)的具体形式为：与仪)小£陶毛s(x)-< Sz&IxeXxJ其中Sj(x)在%上是三次多项式Si(x) = atx3 +瓦/ + ctx + di由插情条件SGD = %, i=0,1,2, ,!),得n+1个条件。边界条件一：ST%) =y(/，s<如)=%/边界条件二：S(%) = %,S(/) = yn边界条件三:假定函数y = f(x)是以b-a为周期的周期函数，这时要求S(x) 也是周期函数，即(SOo + 0) = S(xn - 0) S'(xo

4、 + 0) = Sxn -0)(s (勺 +0) = S"Gn + 0)针对三种边界条件的求解方法的不同，可以分为三转角算法和三弯矩 算法，本实验将介绍和学习三转角算法。三弯矩算法：在每个子区间因为S(x)是三次多项式，所以S(x)是一次多 项式，假设节点处S"Cq) = Mifi = 0,1,2, ,n 则在岛上,/、 x - XiX - Xi_rXi XX Xi.Si(x) =l-M + M1 +1 1 MiXi-1 - Xi一-1hihi对上式积分两次得/、(%i X)3 (X Xj-1)3,、SG) = Mj+ Mt-+ GOi -X)+ 40 勺_1)orii其中

5、c,和4的值为任意常数，由端点值可以求得。显然只要确定出M(i = 0,12，九)就能求出SG),进而得到三次样条函 数S(x)。关于Mj,我们利用样条函数在节点处一阶导数连续来确定。3、算法实例已知函数y = f (x)的如下表的数据 101231245%1312求满足自然边界条件S(x0) = S(m) = 0的”3)的近似值解：根据题意可知，函数在区间两端的二阶导数，因此可知应该用三 弯矩算法来实现，具体程序如下：ttinclude "stdio. h"ttdefine N 4 void main () int i, k;float X,s,yO,yn:float a

6、NN+l, hN, uN, vN, gN, mN, pN, qN, wN;printf ("please input X:");X 为未知数的大小scanf (*%f*,4X);输入x的大小输入y的大小计算步长printf(*please input x:");for(i=0;i<N;i+)scanf (*%f&ai 0);printf(*please input y:");for(i=0:i<N:i+)scanf1);for(i=l;i<N;i+)hi=ai 0-ai-l 0; for(i=l;i<N;i+) vi=hi

7、+l/(hi+hi+l);ui = l-vi;gi=3*ui*(ai+ll)/hi+l+3*vi*(aill)/hi;)printf ("t已知边界条件ln);printf(*t已知边界条件2n");printf ("请选择边界条件序号：");scanf (*%d*,&k);if(k=l)printf ("请输入y0和yn的一阶导：”)；输入边界条件scanf ”, &m0, &mN-l);p0=0;用追赶法求解mNq0=0;gN-2 =gN-2 -u N-2 *mN-l;for (i=l; i<N; i+)(wi

8、=2-ui*pi-l;pi=vi/wi;qi = (gi-ui*qi-l)/wi;)mN-2=qN-2;for (i=N-3; i>0; i一)mi=qi-pi*mi+l;printf ("输出 mi的值:n");for(i=0;i<N;i+)计算最printf (*%fn*, mi); for(i=l;i<N;i+)终结果if(X>ai-l 0&&X<ai 0)s=(hi+2*(X-ai-l 0)*(X-ai 0)*(X-ai 0)*ai-l 1/(hi *hi*hi) + (hi-2*(X-ai 0)*(X-ai-l 0)*

9、(X-ai-l 0)*a i I/(hi*hi*hi) + (X-ai-l 0)*(X-ai 0)*(X-ai 0)*m i-1/ (hi*hi) + (Xai1 0)*(Xai1 0)*(Xai O)*mi/ (hi*hi);printf ("s (%f) =%fn*, X, s);if(k=2)printf ("请输入yO和yn的二阶导：)；输入边界条件二scanf ", &y0, &yn);g0=3*(al l-a0 l)/hl-hl*yO/2;gN-l=3*(aN-l l-aN-2 l)/hN-l+hN-l*yn/2;q0=g0;u0=l;

10、vN-l=l;w0=2;for (i=l; i<N; i+)(wi=2-vi*ui-l/wi-l; qi=gi-vi*qi-l/wi-l;)mN-l=q N-l/wN-l;for (i=N-2; i>=0; i一)mi = (qi-ui*mi+l)/wi;printf ("输出 mi的值:n");for(i=0;i<N;i+)printf mi);for(i=l;i<N;i+)if (X>=ai-1 0ft&X<=ai 0)s=(hi+2*(X-ai-l O)*(X-ai O)*(X-ai O)*ai-1 l/(hi *hi *hi) + (hi-2*(X-ai0)*(X-ai-l0)*(X-ai-l0)*a i l/(hi*hi*hi)4-(X-ai-l 0)*(X-ai 0)*(X-ai 0)*m i-l/(hi*hi) + (X-ai-l 0)*(X-ai-l 0)*(X-ai 0)*mi/ (hi*hi);printf (*s (%f) =%fn*, X, s);)运行结果:please input X：3please input x： 1 2 4 5please input y：1 342(1)己知边界条件1(2)已知边界条件2 请选择边界条件序号：2 请输入W和yn的二阶导：6 0 输出mi的值：2.