牛顿插值法C语言程序
牛顿插值C语言程序
// xg.cpp : 定义控制台应用程序的入口点。
#include"stdafx.h"#include<iostream>#include"math.h"#define N 4using namespace std;void main(){ void lin(double x[],double y[],double t,int n);void newton(double a[],double b[],double t,double h,int n);double t,h,d; int i,j,n,k;double x[N]={0.4,0.5,0.6,0.7};double y[N]={0.38942,0.47943,0.56464,0.64422};double a[N],b[N];h=0.1; //h为等距节点宽度，t为插值点t=0.57891;if(!h) cout<<"此函数为常数"<<endl;else{if(t<x[0]||t>x[N-1])cout<<t<<"不在["<<x[0]<<","<<x[N-1]<<"]范围内"<<endl;else{cout<<"利用分段线性插值："<<endl;lin(x,y,t,N);cout<<"利用等距节点牛顿插值："<<endl;newton(x,y,t,h,N);}}}void lin(double x[],double y[],double t,int n){int i;double w,e,c; //w为逼近值e为余项for(i=0;i<n;i++){if(t==x[i]){cout<<"sin("<<t<<")近似值为"<<y[i]<<endl<<"估计误差为"<<endl;}else if (t>x[i]){ w=y[i]*(t-x[i+1])/(x[i]-x[i+1])+y[i+1]*(t-x[i])/(x[i+1 ]-x[i]);c=fabs(cos(x[i])) > fabs(cos(x[i+1])) ?fabs(cos(x[i])) : fabs(cos(x[i+1]));e=c*(t-x[i])*(t-x[i+1]);cout<<"sin("<<t<<")近似值为"<<w<<endl<<"估计误差为"<<e<<endl;break;}}}void newton(double a[],double b[],double t,double h,int n) {int i,j,c;double d=0.0;double r=1.0,k=0.0,m=1.0,s=0.0; //d为所求近似值，r为差值余项double f(double b[],int n,int a);for(i=0;i<N;i++)if (t<a[i])break;d=b[i-1]; //b[i-1]为牛顿向前插值的第一个点函数值k=(t-a[i-1])/h; //k为（插值点—相邻前一点）/等距hc=j=i-1;for(i=j,j=0;i<N-c;i++,j++){m*=(k-j)/(j+1);if((i+1)<N)s=m * f(b,i,c);d+=s;r*=h*(k-j)/(j+1);}r*=h*(k-j)/(j+1)*fabs(cos(a[c]));cout<<"sin("<<t<<")近似值为"<<d<<endl<<"估计误差为"<<r<<endl;}double f(double b[],int n,int a) //求n阶差分，用p【a】或q 【a】存储n阶差分值{double p[N-1];double q[N-1]; //定义两个数组存储差分int i,j;for(i=0;i<=n;i++){p[i]=b[i+1]-b[i];}for(i=1;i<=n;i++){for(j=0;j<=(N-i-1);j++){if((i%2)==1)q[j]=p[j+1]-p[j];elsep[j]=q[j+1]-q[j];}}if (n%2) return p[a];else return q[a];}/*依次存储n+1阶差分值图示说明： n= 1 2 3 4y[0] p[0] q[0] p[0] q[0]y[1] p[1] q[1] p[1]y[2] p[2] q[2]y[3] p[3]y[4]*/。
c++实现牛顿插值法实验报告概论
数值实验用Newton商差公式进行插值姓名：陈辉学号：13349006院系：数据科学与计算机学院专业：计算机科学与技术班级：计科一班日期：2015-10-11指导老师：纪庆革目录一、实验目的 (3)二、实验题目 (3)三、实验原理与基础理论 (3)四、实验内容 (6)五、实验结果 (10)六、心得体会 (14)七、参考资料 (14)八、附录（源代码） (14)一、实验目的编写一个程序，用牛顿差商公式进行插值。
二、实验题目编写一个程序，用牛顿差商公式进行插值。
三、实验原理与基础理论牛顿插值公式为：N n(x)=f(x0)+f[x0,x1](x−x0)+⋯+f[x0,…,x n](x−x0)…(x−x n−1)其中，f[x0,x1]=f(x1)−f(x0)x1−x0f[x0,…,x k]=f[x1,…x k]−f[x0,…,x k]k0我们将从键盘读入n阶牛顿插值的n+1个节点(x i,f i),i=0,1,…,n，以此得出牛顿插值多项式。
有了节点，我们只需要求f[x0,…,x k]即可。
我们记f[x m,x m−1,…,x m−k]为t[m][k]，则t[m][k]在差商表表的位置为(m, k)：由f[x0,…,x k]的公式可得t[m][k] = (t[m][k-1] - t[m-1][k-1]) / (x[i] – x[i-j])，直观上来看，表中的某格的差商值可由其左边和左上边的值求得，从而牛顿插值多项式变为N(x) = t[0][0] + t[1][1](x-x[0]) + … + t[n][n](x-x[0])…(x-x[n-1])做到这一步，我们已经可以通过上面的数据计算对于给出的x，f(x)的近似值。
为了更规范地表示，下面我把N(x)变换成标准的降幂多项式N(x) = a[n]x^n + a[n-1]x^(n-1) + … + a[2]x^2 + a[1]x + a[0]。
为了便于运算，不妨先把x-x[i]中的减号换成加号（在最后令y[i]=-x[i]，在令x[i]=y[i]，仍可以得到原本的结果），那么有：(x+x[0])(x+x[1])…(x+x[m−1])m−1=x m+x m−1∑x[i]i=0+x m−2∑x[i[0]]x[i[1]]+⋯0≤i[0]≠i[1]≤m−1+x0∑x[i[0]]x[i[1]]…x[i[m−1]]0≤i[0]≠⋯≠i[m−1]≤m−1为了便于表示，把∑x[i[0]]x[i[1]]…x[i[k−1]]0≤i[0]≠⋯≠i[k]≤m−1记为∑x[k]m那么(x+x[0])(x+x[1])…(x+x[m−1])=x m+x m−1∑x[1]+⋯+x0∑x[m]mm只要把N(x)中的每一项展开然后x次数相同的系数相加就可以得到数组a。
实验一牛顿插值法
实验一 牛顿K 次插值多项式一、实验目的：1、掌握牛顿插值法的基本思路和步骤。
2、 培养编程与上机调试能力。
二、牛顿插值法基本思路与计算步骤：给定插值点序列（))(,i i x f x ,,,1,0,n i =。
构造牛顿插值多项式)(u N n 。
输入要计算的函数点,x 并计算)(x N n 的值，利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面)(x N n 的各项系数恰好又是各阶均差，而各阶均差可用均差公式来计算。
为 的 一阶均差。
为的 k 阶均差。
均差表：1． 输入n 值及（))(,i i x f x ,,,1,0,n i =；要计算的函数点x 。
2． 对给定的,x 由[][][]00010101201101()()(),()(),,()()(),,n n n N x f x x x f x x x x x x f x x x x x x x x x f x x x -=+-+--++---计算()n N x 的值。
3．输出()n N x 。
三：程序流程图：四：程序清单：function[c, d]=newpoly(x, y)%牛顿插值的MA TLAB实现%这里x为n个节点的横坐标所组成的向量，y为纵坐标所组成的向量。
%c为所求的牛顿插值多项式的系数构成的向量。
n=length(x);%取x的个数。
d=zeros(n, n);%构造nXn的空数组。
d(: , 1)=y';f or j=2 : nfor k=j : nd(k, j)=(d(k, j-1) - d(k-1, j-1)) / (x(k)-x(k-j+1));ende ndc =d(n, n);for k=(n-1) : - 1 : 1c =conv(c, poly(x(k)));% conv求积，poly(x)将该多项式的系数赋给向量。
m=length(c);c(m)=c(m)+d(k, k);end五、测试数据与结果：测试数据：（第三章习题第三题第2题）01234Y0=-0.916291, y1=-0.693147, y2=-0.510826, y3=-0.357765, y4=-0.223144 建立一个主程序np.mclcclearnewpoly([0.4,0.5,0.6,0.7,0.8],[ -0.916291, -0.693147, -0.510826, -0.357765, -0.223144]) 计算结果如下：ans =-0.3096 2.6083 -5.4861 5.6921 -2.4744由此看出所求的牛顿多项式为：P(x)= -0.3096x4+2.6083x3-5.4861x2+5.6921x-2.4744P(0.53)= -0.6347。
用C语言实现牛顿向前插值计算
用C语言实现牛顿向前插值计算，程序代码如下：#include "stdlib.h"#include "stdio.h"#include "conio.h"#include "string.h"#include "math.h"#define N 100typedef struct{float x;float y;}POINT;float CreTable(int n,POINT Table[N],float y[N][N]){int i,j,count=0;for(i=0;i<n;i++){printf("Input %dth point(x%d,y%d):",i+1,i+1,i+1);scanf("%f,%f",&Table[i].x,&Table[i].y);}for(i=0;i<n;i++)y[i][0]=Table[i].y;for(j=1;j<n;j++,count++){i=count;for(;i<n-1;i++)y[i+1][j]=y[i+1][j-1]-y[i][j-1];}}float ShowTable(int n,POINT Table[N],float y[N][N]){int i,j;printf("\nForward difference table:\n");printf(" x y\n");for(i=0;i<n;i++){printf("%10.5f",Table[i].x);for(j=0;j<=i;j++){printf("%10.5f",y[i][j]);if(j==i) printf("\n");}}}float Calculus(int n,POINT Table[N],float y[N][N]) {int i,j,k,count[N]={0};float tmp1,tmp2,tmp3,sum,a[N],h[N],x[N],t[N];printf("\nNumber of x:");scanf("%d",&k);for(j=0;j<k;j++){printf("Input x0[%d]:",j+1);scanf("%f",&a[j]);printf("Input x[%d]:",j+1);scanf("%f",&x[j]);printf("Input h[%d]:",j+1);scanf("%f",&h[j]);t[j]=(x[j]-a[j])/h[j];printf("t[%d]=%.5f\n",j+1,t[j]);for(i=0;i<n;i++)if(Table[i].x==a[j]){count[j]=i;break;}}printf("\nThe results:");for(j=0;j<k;j++){tmp1=1;tmp2=1;sum=y[count[j]][0];for(i=1;;i++){tmp1=tmp1*(t[j]-i+1);tmp2=tmp2*i;tmp3=tmp1*y[count[j]+i][i]/tmp2;sum=sum+tmp3;if(count[j]+i==n-1) break;}printf("\n N(%f)= %f",x[j],sum);}}int main(){int n;float y[N][N],x[N];POINT Table[N];printf("Number of points n=");scanf("%d",&n);CreTable(n,Table,y);ShowTable(n,Table,y);Calculus(n,Table,y);getch();}。
牛顿插值法的C语言编程
Newton 插值Newton 插值函数 Newton 插值函数是用差商作为系数，对于01,,,n x x x …这1n +个点，其一般形式为：00100120101011()[][,]()[,,]()()[,,,]()()()n n n N x f x f x x x x f x x x x x x x f x x x x x x x x x −=+−+−−++−−−…………对于011,,,n x x x −…这n 个点，100100120101012()[][,]()[,,]()()[,,,]()()()n n n N x f x f x x x x f x x x x x x x f x x x x x x x x x −−=+−+−−++−−−…………差商的定义 若已知函数()f x 在点(0,1,2,,)i x i n =⋅⋅⋅处的函数值()i f x 。
则称：00[]()f x f x =为函数()f x 在点0x 的0阶差商；100110[][][,]f x f x f x x x x −=−为函数()f x 关于01,x x 的1阶差商；120101220[,][,][,,]f x x f x x f x x x x x −=−为函数()f x 过点012,,x x x 的2阶差商；依此类推，一般地称121012101210[,,,,][,,,,][,,,,,]k k k k k k k f x x x x f x x x x f x x x x x x x −−−−⋅⋅⋅−⋅⋅⋅⋅⋅⋅=−为函数()f x 关于01,,,k x x x ⋅⋅⋅的k 阶差商。
表1 差商表i x ()i f x 1阶差商 2阶差商3阶差商4阶差商0x 1x 2x 3x 4x……0()f x 1()f x 2()f x 3()f x 4()f x ……01[,]f x x12[,]f x x 23[,]f x x 34[,]f x x……012[,,]f x x x 123[,,]f x x x 234[,,]f x x x ……0123[,,,]f x x x x 1234[,,,]f x x x x ……01234[,,,,]f x x x x x……根据Newton 插值函数编写的C 语言编程根据Newton 插值函数并对照上面的差商表，可编写出Newton 插值法的C 语言程序如下： #include<stdio.h> #include<iostream.h> #include <stdlib.h>double NewtonInterpolation(double *x,double *y,int n,double xx,double *pyy) {double *f=(double *)malloc(n*sizeof(double));int i,k;for(i=1;i<=n-1;i++)f[i]=y[i];for(k=1;k<=n-1;k++)for(i=k;i<=n-1;i++){f[i]=(y[i]-y[i-1])/(x[i]-x[i-k]);if(i==n-1)for(i=k;i<n;i++)y[i]=f[i];}*pyy=y[n-1];for(i=n-2;i>=0;i--)*pyy=(*pyy)*(xx-x[i])+y[i];free(f);return 0;}void main(){int n=5;double x[5]={1.0,2.7,3.2,4.8,5.6},y[5]={14.2,17.8,22.0,38.3,51.7},xx=3,yy; NewtonInterpolation(x,y,n,xx,&yy);printf("%lf\n",yy);}。
Newton插值和三次样条插值的程序代码
(){}21()(11),5,10,20:12521()1,(0,1,2,,)()2,(0,1,2,,)()()235,20:1100(),i i i i n n k k k Newton f x x n x f x x i i n f x n x y i n Newton N x S x n x k y f x N =-≤≤=+=-+====-+=插值多项式和三次样条插值多项式。
已知对作、计算函数在点处的值；、求插值数据点的插值多项式和三次样条插值多项式；、对计算和相应的函数值()() (1,2,,99)4:()max()()max ()n k n k n k n k n k n k k k x S x k E N y N x E S y S x ==-=-和；、计算，；解释你所得到的结果。
Newton 插值的程序代码：%求Newton 插值多项式function new=newton(m,y)n=length(m);b=1;new=[zeros(1,n-1),y(1)];for i=2:n y(i:n)=(y(i:n)-y(i-1:n-1))./(m(i:n)-m(1:n-i+1)); %求差商b=conv(b,[1,-m(i-1)]);a=[zeros(1,n-length(b)),b];new=new+y(i)*a;end三次样条插值的程序代码：%求三次样条插值多项式function s=myspline(m,y)syms x;n=length(m);M=y;M(2:n)=(M(2:n)-M(1:n-1))./(m(2:n)-m(1:n-1));M(3:n)=(M(3:n)-M(2:n-1))./(m(3:n)-m(1:n-2));M=[0,M(3:n),0];h=m(2:n)-m(1:n-1);c=h(2:n-1)./(h(2:n-1)+h(1:n-2));a=1-c;b=2+zeros(1,n);M=6*M;c=[0,c];a=[a,0];M=chase(a,b,c,M); %用追赶法解方程组for i=1:n-1 %分区间计算插值多项式s(i)=M(i)*(m(i+1)-x)^3+M(i+1)*(x-m(i))^3+(6*y(i)-M(i)*h(i)^2)*(m(i+1)-x)+(6*y(i+1)-M(i+1)*h(i)^2 )*(x-m(i));s(i)=s(i)/(6*h(i));end主程序为：clear;clc;n=5; %指定n的值syms x %符号数xf=1./(1+25*x^2); %定义函数f(x)i=0:n;m=-1+2*i/n;y=subs(f,m); %计算函数值f(x)disp(['When n=',num2str(n),',the values of f(x) at the points of x(i) are:']);disp(y);%求Newton插值多项式new=newton(m,y);new=round(10000*new)/10000;nout=poly2sym(new);nout=vpa(nout,4); %整理多项式并输出disp(['When n=',num2str(n),',N(x)=']);disp(nout);%求三次样条插值多项式s=myspline(m,y);for i=1:n %分区间计算插值多项式a=sym2poly(s(i));b=poly2sym(a);b=vpa(b,4); %整理多项式并输出disp(['In the range of [',num2str(m(i)),',',num2str(m(i+1)),'],S(x)=']);disp(b);endk=0:100;xk=-1+0.02*k; %计算xk的值ll1=subs(f,xk);ll2=polyval(new,xk); %计算f(xk),N(xk)for i=1:length(xk)r=find(m(2:n+1)>=xk(i),1,'first');ll3(i)=subs(s(r),xk(i)); %计算S(xk)enden=max(abs(ll1-ll2));es=max(abs(ll1-ll3)); %计算E(Nn),E(Sn),并输出disp(['E(Nn)=',num2str(en)]);disp(['E(Sn)=',num2str(es)]);plot(xk,ll1,xk,ll2,'r',xk,ll3,'g');axis([-1,1,-0.1,1.1]);title('Interpolation');legend('Original function','Newton Interpolation','Spline')grid; %将原函数，Newton插值函数，三次样条插值函数画成图形。
拉格朗日和牛顿插值法的C 方法实现(数值分析上机实验)
数值分析上机实验实验一一．上机题目：已知： 4 =2,9 =3,16 =4分别用二次Lagrange和Newton插值法求7 的近似值。
二．解题方法：1．lagrange方法：设x0=4，y0=2，x1=9，y1=3，x2=16，y2=4代入方程:(x1-X)(x2-X)/(x1-x0)(x2-x0)*y0+(x0-X)(x2-X)/(x0-x1)(x2-x1)*y1+(x1-X)(x0-X)/(x1-x2)(x0-x2)*y2令X=7代入方程得 Y=2.628572．Newton方法：设x0=4，y0=2，x1=9，y1=3，x2=16，y2=4建表4 29 3 0.216 4 0.14286 -0.00476f(x)=f(x0)+f[x0,x1](X-x0)+f[x0,x1,x2](X-x0)(X-x1)(X-x2)令X=7代入方程得Y=2.62857三．算法公式步骤：grange方法：通过公式写出算法并得出最后的值Y：for(b=0;b<m;b++)//完成公式f(Xn)外层嵌套循环f[b]=i//{double l=1;//保证每次跳出内层循环将L置1 不会将第一项的值带入下一项//for(a=0;a<m;a++)//完成公式f(Xn)内层嵌套循环f[a]=j//{if(a!=b)//完成定义i=1,i!=j//l=(f[a]-F)/(f[a]-f[b])*l;//完成(j-m)/(j-i)//la=l*g[b];//完成公式的F(X0)=f(X0)*Y0并累乘输出结果// }Y=la+Y;//累加x0y0+x1y1+...得最后结果//}2．Newton方法：先建表，通过二维数组的思想建表for(l=2;l<m+2;l++)//外层循环控制y阶数//{for(k=1;k<m+1;k++)//内层循环控制x个数//{a[k][l]=(a[k][l-1]-a[k-1][l-1])/(a[k][0]-a[k-l+1][0]);//完成f(x0,x1,...,xn)并存表//}}填表。
牛顿插值函数C语言程序实现
牛顿插值函数C语言程序实现牛顿插值的关键在于差商表的计算，差商表第一行是y值，为了配合计算，在该矩阵上方配上节点x0、x1、x2……xnf[x0,x1]=[f(x1)−f(x0)]/(x1−x0)f[x0,x1]=[f(x1)−f(x0)]/(x1−x0)f[x0,x1,x2]=[f[x1,x2]−f[x0,x1]]/(x2−x0)f[x0,x1,x2]=[f[x1,x2]−f[x0,x1]]/(x2−x0)……所以只要计算矩阵内上三角值即可。
#include <stdio.h>#include <stdlib.h>int main(){float table(int n,float a1[10],float a2[10],float a3[10][10]);float newton(int n,float a4[10][10],float a5[10]);float arrX[10],arrY[10],arrL[10][10];int num,i;printf("请输入插值节点的个数（个数应小于10）：");scanf("%d",&num);printf("请输入各个插值节点的值：\n");for(i=0; i<num; i++){printf("请输入X%d值：",i+1);scanf("%f",&arrX[i]);printf("请输入Y%d值：",i+1);scanf("%f",&arrY[i]);}table(num,arrX,arrY,arrL);newton(num,arrL,arrX);return 0;}float table(int n,float a1[10],float a2[10],float a3[10][10]){int i,j;for(i=0; i<n; i++){a3[0][i]=a2[i];//第一行初始化为y值}for(i=0; i<n; i++)for(j=n-1; j>i; j--)//从一行最后往前循环，到i=j为止，即上三角全部计算赋值{a3[i+1][j]=(a3[i][j]-a3[i][j-1])/(a1[j]-a1[j-1*(i+1)]);//差商表计算，最后一项arrX[j-1*(i+1)]脚标是计算步长}for(i=1; i<n; i++)//下三角未计算赋值，系统会随机分配值给下三角的每一位，故赋值0会使差商表输出更整齐for(j=0; j<i; j++){a3[i][j]=0.0;}printf("差商表为：\n");printf("----------------------------------------------------------\n");for(i=0; i<n; i++){for(j=0; j<n; j++){if(j%n==0)//num个数一行输出printf("\n");printf("%f\t",a3[i][j]);}}printf("\n");printf("----------------------------------------------------------\n");return 0;}float newton(int n,float a4[10][10],float a5[10]){int i;float x,y,t1=1.0;while(1){printf("请输入要插入节点的X值：");scanf("%f",&x);y=a4[0][0];for(i=1; i<n; i++){t1=t1*(x-a5[i-1]);y=y+a4[i][i]*t1;//差商表对角线值依次乘以(x-x0)(x-x1)……}printf("插值结果为：%f",y);printf("\n");}return 0;}运行结果如下：。