三次樣條函數(shù)三彎矩算法

1、摘要求函數(shù)在給定區(qū)間上的定積分，在微積分學(xué)中已給出了許多計(jì)算方法, 但是，在實(shí)際問題計(jì)算中，往往僅給出函數(shù)在一些離散點(diǎn)的值，它的解析 表達(dá)式?jīng)]有明顯的給出，或者，雖然給出解析表達(dá)式，但卻很難求得其原 函數(shù)，這時(shí)我們可以通過數(shù)值方法求出函數(shù)積分的近似值。在用近似值代替真實(shí)值時(shí)，遇到的問題就是近似值的代數(shù)精度是否足 夠。當(dāng)代數(shù)精度不足夠時(shí)，很顯然提高插值函數(shù)的次數(shù)是一種方法，但是 考慮到數(shù)值計(jì)算的穩(wěn)定性，當(dāng)次數(shù)過高時(shí)，會出現(xiàn)龍格現(xiàn)象，用增大n的 方法來提高數(shù)值積代數(shù)精度是不可取的。正如我們所知道的分段線性插值，逼近程度好，但光滑性差。分段三 次Hermite插值，逼近程度好，光滑性也有所提高，但是

2、需要增加更多的 條件，不太實(shí)用。因此，我們將介紹一種結(jié)合二者的優(yōu)點(diǎn)的插值方法一一 三次樣條插值。本實(shí)驗(yàn)將介紹三次樣條插值的三彎矩算法。關(guān)鍵詞：龍格現(xiàn)象三彎矩算法代數(shù)精度分段三次Hermite插值1、實(shí)驗(yàn)?zāi)康?)通過本次實(shí)驗(yàn)體會并學(xué)習(xí)三次樣條插值的優(yōu)點(diǎn)。2)通過對三次樣條插值進(jìn)行編程實(shí)現(xiàn)，提高自己的編程能力。3)用實(shí)驗(yàn)報(bào)告的形式展現(xiàn)，提高自己在寫論文方面的能力。2、算法流程如果已知函數(shù)y= f(x)在節(jié)點(diǎn)a = %o <<V = b,yt = 0,1,2,，幾處的函數(shù)值和導(dǎo)數(shù)值：y = f("),i = 0,LAm如果S(x)滿足條件:1) S(x)是一個(gè)分段的三次多項(xiàng)式且

3、*(x) = v ；2) S(x)在a,b具有二階連續(xù)導(dǎo)數(shù)。則稱S(x)是三次樣條插值函數(shù)。S(x)的具體形式為：與儀)小£陶毛s(x)-< Sz&IxeXxJ其中Sj(x)在%上是三次多項(xiàng)式Si(x) = atx3 +瓦/ + ctx + di由插情條件SGD = %, i=0,1,2, ,!),得n+1個(gè)條件。邊界條件一：ST%) =y(/，s<如)=%/邊界條件二：S(%) = %,S(/) = yn邊界條件三:假定函數(shù)y = f(x)是以b-a為周期的周期函數(shù)，這時(shí)要求S(x) 也是周期函數(shù)，即(SOo + 0) = S(xn - 0) S'(xo

4、 + 0) = Sxn -0)(s (勺 +0) = S"Gn + 0)針對三種邊界條件的求解方法的不同，可以分為三轉(zhuǎn)角算法和三彎矩 算法，本實(shí)驗(yàn)將介紹和學(xué)習(xí)三轉(zhuǎn)角算法。三彎矩算法：在每個(gè)子區(qū)間因?yàn)镾(x)是三次多項(xiàng)式，所以S(x)是一次多 項(xiàng)式，假設(shè)節(jié)點(diǎn)處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的值為任意常數(shù)，由端點(diǎn)值可以求得。顯然只要確定出M(i = 0,12，九)就能求出SG),進(jìn)而得到三次樣條函 數(shù)S(x)。關(guān)于Mj,我們利用樣條函數(shù)在節(jié)點(diǎn)處一階導(dǎo)數(shù)連續(xù)來確定。3、算法實(shí)例已知函數(shù)y = f (x)的如下表的數(shù)據(jù) 101231245%1312求滿足自然邊界條件S(x0) = S(m) = 0的”3)的近似值解：根據(jù)題意可知，函數(shù)在區(qū)間兩端的二階導(dǎo)數(shù)，因此可知應(yīng)該用三 彎矩算法來實(shí)現(xiàn)，具體程序如下：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 為未知數(shù)的大小scanf (*%f*,4X);輸入x的大小輸入y的大小計(jì)算步長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的一階導(dǎo)：”)；輸入邊界條件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+)計(jì)算最printf (*%fn*, mi); for(i=l;i<N;i+)終結(jié)果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的二階導(dǎo)：)；輸入邊界條件二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);)運(yùn)行結(jié)果:please input X：3please input x： 1 2 4 5please input y：1 342(1)己知邊界條件1(2)已知邊界條件2 請選擇邊界條件序號：2 請輸入W和yn的二階導(dǎo)：6 0 輸出mi的值：2.