/*???三次樣條插值計算算法???*/???
??#include???“math.h“???
??#include???“stdio.h“???
??#include???“stdlib.h“???
??
??/*???
??N：已知節(jié)點數(shù)N+1???
??R：欲求插值點數(shù)R+1???
??x，y為給定函數(shù)f（x）的節(jié)點值{x（i）}???（x（i）??P0=f（x0）的二階導(dǎo)數(shù)；Pn=f（xn）的二階導(dǎo)數(shù)???
??u:存插值點{u（i）}???????0<=i<=R???
??求得的結(jié)果s（ui）放入s[R+1]???????0<=i<=R???
??返回0表示成功，1表示失敗???
??*/
??
??int???SPL（int???Nint???Rdouble???x[]double???y[]double???P0double???Pndouble???u[]double???s[]）???
??{???
??/*聲明局部變量*/???
??double???*h;???/*存放步長:{hi}???????0<=i<=N-1?????*/???
??double???*a;???/*存放系數(shù)矩陣{ai}???????1<=i<=N???；???分量0沒有利用?????*/???
??double???*c;???/*先存放系數(shù)矩陣{ci}?????后存放{Bi}?????0<=i<=N-1?????*/???
??double???*g;???/*先存放方程組右端項{gi}?????后存放求解中間結(jié)果{yi}?????0<=i<=N?????*/???
??double???*af;???/*存放系數(shù)矩陣{a（f）i}???????1<=i<=N???；?????*/???
??double???*ba;???/*存放中間結(jié)果???0<=i<=N-1*/???
??double???*m;???/*存放方程組的解{m（i）}???????0<=i<=N???；?????*/???
????
??int???ik;???
??double???p1p2p3p4;???????
????
??/*分配空間*/???
??if（!（h=（double*）malloc（N*sizeof（double））））???exit（1）;???
??if（!（a=（double*）malloc（（N+1）*sizeof（double））））???exit（1）;???
??if（!（c=（double*）malloc（N*sizeof（double））））???exit（1）;???
??if（!（g=（double*）malloc（（N+1）*sizeof（double））））???exit（1）;???
??if（!（af=（double*）malloc（（N+1）*sizeof（double））））???exit（1）;???
??if（!（ba=（double*）malloc（（N）*sizeof（double））））???exit（1）;???
??if（!（m=（double*）malloc（（N+1）*sizeof（double））））???exit（1）;???
????
??/*第一步：計算方程組的系數(shù)*/???
??for（k=0;k??h[k]=x[k+1]-x[k];???
??for（k=1;k??a[k]=h[k]/（h[k]+h[k-1]）;???
??for（k=1;k??c[k]=1-a[k];???
??for（k=1;k??g[k]=3*（c[k]*（y[k+1]-y[k]）/h[k]+a[k]*（y[k]-y[k-1]）/h[k-1]）;???
??c[0]=a[N]=1;???
??g[0]=3*（y[1]-y[0]）/h[0]-P0*h[0]/2;???
??g[N]=3*（y[N]-y[N-1]）/h[N-1]+Pn*h[N-1]/2;???
????
??/*第二步：用追趕法解方程組求{m（i）}?????*/???
??ba[0]=c[0]/2;???
??g[0]=g[0]/2;???
??for（i=1;i??{???
??af[i]=2-a[i]*ba[i-1];???
??g[i]=（g[i]-a[i]*g[i-1]）/af[i];???
??ba[i]=c[i]/af[i];???
??}???
??af[N]=2-a[N]*ba[N-1];???
??g[N]=（g[N]-a[N]*g[N-1]）/af[N];???
????
??m[N]=g[N];???????????????????/*P110???公式：6.32*/???
??for（i=N-1;i>=0;i--）???
??m[i]=g[i]-ba[i]*m[i+1];???
????
??/*第三步：求值*/???
??for（i=0;i<=R;i++）???
??{???
??/*判斷u（i）屬于哪一個子區(qū)間，即確定k???*/???
??if（u[i]x[N]）???????
??{???
??/*釋放空間*/???
??free（h）;???
??free（a）;???
??free（c）;???
??free（g）;???
??free（af）;???
??free（ba）;???
??free（m）;???
??return???1;???
??}???
??k=0;???
??while（u[i]>x[k+1]）???
??k++;???
????
????
??p1=（h[k]+2*（u[i]-x[k]）*pow（（u[i]-x[k+1]）2）*y[k]）/pow（h[k]3）;???
??p2=（h[k]-2*（u[i]-x[k+1]）*pow（（u[i]-x[k]）2）*y[k+1]）/pow（h[k]3）;???
??p3=（u[i]-x[k]）*pow（（u[i]-x[k+1]）2）*m[k]/pow（h[k]2）;???
??p4=（u[i]-x[k+1]）*pow（（u[i]-x[k]）2）*m[k+1]/pow（h[k]2）;???
??s[i]=p1+p2+p3+p4;???
??}???
????
??/*釋放空間*/???
??free（h）;???
??free（a）;???
??free（c）;???
??free（g）;???
??free（af）;???
??free（ba）;???
??free（m）;???
????
??return???0;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????4792??2017-10-17?23:31??spline\spline.c

?????目錄??????????0??2017-10-17?23:32??spline

-----------?---------??----------?-----??----

?????????????????4792????????????????????2