牛顿插值法的C语言编程

Newton插值的C++实现Newton（⽜顿）插值法具有递推性，这决定其性能要好于Lagrange（拉格朗⽇）插值法。

其重点在于差商（Divided Difference）表的求解。

步骤1. 求解差商表，这⾥采⽤⾮递归法（看着挺复杂挺乱，这⾥就要⾃⼰动笔推⼀推了，闲了补上其思路），这样，其返回的数组（指针）就是差商表了，/** 根据插值节点及其函数值获得差商表* 根据公式⾮递归地求解差商表* x: 插值节点数组* y: 插值节点处的函数值数组* lenX: 插值节点的个数* return: double类型数组*/double * getDividedDifferenceTable(double x[], double y[], int lenX) {double *result = new double[lenX*(lenX - 1) / 2];for (int i = 0; i < lenX - 1; i++) {result[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);}int step = lenX - 1; // 增加步长int yindex = lenX - 1; // 分⼦的基准值，线性减速递增int xgap = 2; // 分母的间距，每次+1int xstep = 2;while (step >= 1) {for (int i = 0; i < step - 1; i++) {result[yindex + i] = (result[yindex - step + i + 1] - result[yindex - step + i]) / (x[xstep + i] - x[xstep + i - xgap]);}yindex += (step - 1);xstep++;step--;xgap++;}return result;}步骤 2. 取得差商表每⼀列的第⼀个作为基函数的系数/*** 从差商表中获取⼀定阶数的某个差商* dd: 差商表* len: 差商表的长度* rank: 所取差商的阶数* index: rank阶差商的第index个差商* return: 返回double类型所求差商*/double getDividedDifference(double dd[], int len, int rank, int index) {// 根据差商表的长度求解差商的最⾼阶数// 由于差商表是三⾓形的，因此有规律可循，解⽅程即可int rankNum = (int)(0.5 + sqrt(2 * len + 0.25) - 1);printf("%d 阶差商 \n", rank);int pos = 0;int r = 1;// 根据n+(n+1)+...+(n-rank)求得rank阶差商在差商表数组中的索引while (rank > 1){// printf("geting dd's index, waiting...\n");pos += rankNum;rank--;rankNum--;}pos += index; // 然后加上偏移量，意为rank阶差商的第index个差商return dd[pos];}步骤3. Newton 插值主流程/** Newton Interpolating ⽜顿插值主要代码* x: 插值节点数组* y: 插值节点处的函数值数组* lenX: 插值节点个数（数组x的长度）* newX: 待求节点的数组* lennx: 待求节点的个数（数组lenX的长度）*/double *newtonInterpolation(double x[], double y[], int lenX, double newX[], int lennx) {// 计算差商表double *dividedDifferences = getDividedDifferenceTable(x, y, lenX);//double *result = new double[lennx];// 求差商表长度int ddlength = 0;for (int i = 1; i < lenX; i++){ddlength += i;}// 打印差商表信息printf("=======================差商表的信息=======================\n");printf("差商个数有：%d, 待插值节点个数：%d\n =======================差商表=======================", ddlength, lennx); int ifnextrow = 0; // 控制打印换⾏for (int i = 0; i < ddlength; i++) {if (ifnextrow == (i)*(i + 1) / 2)printf("\n");printf("%lf, ", dividedDifferences[i]);ifnextrow++;}printf("\n");// 计算Newton Interpolating 插值double ddi = 0;double coef = 1;for (int i = 0; i < lennx; i += 2) {coef = (double)1;result[i] = y[i];printf("=============打印计算过程=============\n");for (int j = 1; j < lenX -1; j++) {coef *= (newX[i] - x[j - 1]);ddi = getDividedDifference(dividedDifferences, ddlength, j, 0);printf("取得差商： %lf\n", ddi);result[i] = result[i] + coef * ddi;printf("计算result[%d] + %lf * %lf", i, coef, ddi);printf("获得结果result %d: %lf\n", i, result[i]);}printf("result[%d] = %lf\n", i, result[i]);// =======================选做题：求误差==========================printf("求解截断误差所需差商:%lf\n", getDividedDifference(dividedDifferences, ddlength, lenX-1, 0));// printf("求解截断误差所需多项式：%lf\n", coef);printf("求解截断误差所需多项式的最后⼀项：%lf - %lf\n", newX[i] , x[lenX - 2]);result[i + 1] = getDividedDifference(dividedDifferences, ddlength, lenX - 1, 0)*coef*(newX[i] - x[lenX - 2]);// ===============================================================}if (dividedDifferences) {delete[] dividedDifferences;}return result;}步骤4. 主函数⽰例int main(){std::cout << "Hello World!\n";printf("Newton插值!\n");double x[] = {0.40, 0.55, 0.65, 0.80, 0.90, 1.05};//double x[] = { 1.5, 1.6, 1.7 };int lenX = sizeof(x) / sizeof(x[0]);double y[] = {0.41075, 0.57815, 0.69675, 0.88811, 1.02652, 1.25382};//double y[] = { 0.99749, 0.99957, 0.99166 };double newx[] = {0.596};//double newx[] = { 1.609 };// 计算⽜顿插值结果和截断误差double *result = newtonInterpolation(x, y, lenX, newx, 1);printf("=====================最终结果=====================\n");printf("求得sin(1.609)的近似值为：%lf\n", *result);printf("截断误差：%.10lf\n", *(result + 1));if (result) {delete[] result;}return0;}。

实验报告：数学与统计学系信息与计算科学专业实验报告一、题目1、上机作业题程序12、上机作业题程序2二、算法1、Lagrange 插值//输入被插值点的数目POINT;int main(){int n;inti,j;POINT points[MAX_N+1];double diff[MAX_N+1];doublex,tmp=0,lagrange=0,tx,ty;printf("\nInput n value:");scanf("%d",&n);if(n>MAX_N){printf("The input n is larger thenMAX_N,please redefine the MAX_N.\n");return 1;}if(n<=0){printf("Please input a number between 1 and %d\n",MAX_N);return 1;}//输入被插值点printf("Now input the (x_i,y_i),i=0,...,%d:\n",n);for(i=0;i<=n;i++)scanf("%lf%lf",&points[i].x,&points[i].y);printf("Now input the x value:"); //输入计算Lagrange插值多项式的x值scanf("%lf",&x);for(i=0;i<=n;i++){diff[i]=0;tx=1;ty=1;for(j=0;j<=n;j++){if(i!=j){tx=tx*(x-points[j].x);ty=ty*(points[i].x-points[j].x);}}diff[i]=tx/ty;}for(i=0;i<=n;i++){tmp=points[i].y*diff[i];printf("%f",tmp);lagrange+=tmp;}printf("lagrange(%f)=%f\n",x,lagrange);return 0;}2、Newton 插值//输入被插值点的数目POINT;int main(){ int n;inti,j;POINT points[MAX_N+1];double diff[MAX_N+1];doublex,tmp,newton=0;printf("\nInput n value: ");scanf("%d",&n);if (n>MAX_N){printf("The input n is larger thenMAX_N,please redefine the MAX_N.\n");return 1;}if(n<=0){printf("Please input a number between 1 and %d\n",MAX_N);// getch(); return 1;}//输入被插值点printf("Now input the (x_i,y_i),i=0,...,%d:\n",n);for (i=0;i<=n;i++)scanf("%lf%lf",&points[i].x,&points[i].y);printf("Now input the x value: ");//输入计算Newton插值多项式的x值scanf("%lf",&x);for (i=0;i<=n;i++)diff[i]=points[i].y;for (i=0;i<n;i++){for (j=n;j>i;j--){diff[j]=(diff[j]-diff[j-1])/(points[j].x-points[j-1-i].x);}//计算f(x_0,…,x_n)的差商}tmp=1;newton=diff[0];for(i=0;i<n;i++){tmp=tmp*(x-points[i].x);newton=newton+tmp*diff[i+1];}printf("newton(%f)=%f\n",x,newton);return 0;}三、C程序1、Lagrange 插值#include <stdio.h>#define MAX_N 20typedefstructtagPOINT{double x;double y;}POINT;int main(){int n;inti,j;POINT points[MAX_N+1];double diff[MAX_N+1];doublex,tmp=0,lagrange=0,tx,ty;printf("\nInput n value:");scanf("%d",&n);if(n>MAX_N){printf("The input n is larger thenMAX_N,please redefine the MAX_N.\n");return 1;}if(n<=0){printf("Please input a number between 1 and %d\n",MAX_N); return 1;}printf("Now input the (x_i,y_i),i=0,...,%d:\n",n);for(i=0;i<=n;i++)scanf("%lf%lf",&points[i].x,&points[i].y);printf("Now input the x value:");scanf("%lf",&x);for(i=0;i<=n;i++){diff[i]=0;tx=1;ty=1;for(j=0;j<=n;j++){if(i!=j){tx=tx*(x-points[j].x);ty=ty*(points[i].x-points[j].x);}}diff[i]=tx/ty;}for(i=0;i<=n;i++){tmp=points[i].y*diff[i];printf("%f",tmp);lagrange+=tmp;}printf("lagrange(%f)=%f\n",x,lagrange);return 0;}2、Newton 插值#include <stdio.h>#define MAX_N 20typedefstructtagPOINT{ double x;double y;} POINT;int main(){ int n;inti,j;POINT points[MAX_N+1];double diff[MAX_N+1];doublex,tmp,newton=0;printf("\nInput n value: ");scanf("%d",&n);if (n>MAX_N){printf("The input n is larger thenMAX_N,please redefine the MAX_N.\n");return 1;}if (n<=0){printf("Please input a number between 1 and %d.\n",MAX_N);return 1;}//输入被插值点（x_i,y_i）printf("Now input the (x_i,y_i),i=0,...,%d:\n",n);for (i=0;i<=n;i++)scanf("%lf%lf",&points[i].x,&points[i].y);printf("Now input the x value: ");scanf("%lf",&x);for (i=0;i<=n;i++)diff[i]=points[i].y;for (i=0;i<n;i++){for (j=n;j>i;j--){diff[j]=(diff[j]-diff[j-1])/(points[j].x-points[j-1-i].x);}}tmp=1;newton=diff[0];for(i=0;i<n;i++){tmp=tmp*(x-points[i].x);newton=newton+tmp*diff[i+1];}printf("newton(%f)=%f\n",x,newton);return 0;}四、运行结果1、Lagrange 插值1910年Larange插值计算得到的人口数：1965年Larange插值计算得到的人口数：2002年Larange插值计算得到的人口数：从插值计算得出的结果1910年的人口数是31872000人，1965年的人口数约为193081511人，2002年的人口数约为26138748，而1910年的实际人口数为91772000人，1960年的实际人口数为179323000人,1970年的人口数为203212000人，所以拉格朗日插值计算得出的结果只有1965年的人口数与实际值相差较近，而1910年和2002年的计算结果都与实际值相差较大，所以插值计算得到的数据准确性并不高。

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);}。

(){}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语言牛顿迭代法牛顿迭代法（Newton-Raphson法）是一种求解方程近似解的方法，它是利用泰勒级数展开函数在某点的值，然后用一阶泰勒展开式的根近似表示函数的零点，因此也被称为牛顿拉弗森法。

它可以高效地解决复杂的非线性方程组，是科学计算领域中最为常用和基础的方法之一。

牛顿迭代法的基本思想是：在第k次迭代时，求出曲线f(x)在点xk的一次导数斜率，以此确定x轴上的一个点xk+1，和该点处曲线的一次切线。

这条切线和x轴交点的横坐标就是极值点的估计值。

这个过程可以迭代多次，直到达到满足一定的误差精度或者迭代次数的要求。

C语言实现牛顿迭代法需要先定义一个函数，这个函数就是需要求解方程的函数。

定义完函数之后，需要实现牛顿迭代公式来求出下一次迭代的估计值，然后不断迭代。

具体实现过程如下：1. 定义函数f(x)，即需要求解方程的函数。

2. 定义函数f_prime(x)，即f(x)的一次导数。

3. 定义变量x和x_next，初始化它们的值。

4. 在循环中，首先计算f(x)和f_prime(x)，然后计算下一个迭代点的估计值x_next = x - f(x) / f_prime(x)。

5. 如果x_next和x的差异满足预设的精度要求，则退出循环。

6. 否则，将x_next的值赋值给x，并重复执行第4步。

C语言实现牛顿迭代法的代码如下：#include <stdio.h>#include <math.h>定义函数f(x)double f(double x) {return x * x - 2;}定义函数f_prime(x)double f_prime(double x) {return 2 * x;}int main() {定义变量double x, x_next, epsilon;int iter;初始化变量x = 1.0;epsilon = 1e-6;iter = 0;迭代求解do {x_next = x - f(x) / f_prime(x);iter++;printf("Iteration %d: x = %lf\n", iter, x_next);x = x_next;} while (fabs(x_next - x) >= epsilon);输出结果printf("Final result: x = %lf\n", x);return 0;}在这个代码中，我们使用了do-while循环来不断执行迭代过程，直到达到预设的精度要求。

数值分析上机实验实验一一．上机题目：已知： 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语言实现牛顿法求方程的根的示例代码：```c#include <stdio.h>#include <math.h>// 要求解的方程double f(double x) {return x*x - 2; // 求解x^2 - 2 = 0的根}// 方程的导数double df(double x) {return 2*x; // 方程的导数为2x}// 牛顿法求根函数double newton(double x0, double epsilon) {double x = x0; // 初始值double delta; // 误差do {double fx = f(x);double dfx = df(x);delta = fx / dfx; // 迭代公式x = x - delta; // 迭代计算新的x} while (fabs(delta) > epsilon); // 判断误差是否小于给定的精度return x;}int main() {double x0 = 1; // 初始值double epsilon = 0.0001; // 精度double root = newton(x0, epsilon);printf("方程的根为：%lf\n", root);return 0;}```以上代码中，`f`函数表示要求解的方程，`df`函数表示方程的导数，`newton`函数是用牛顿法求方程的根的函数，`main`函数是主函数，用于测试求解结果。

在代码中，初始值`x0`和精度`epsilon`可以根据实际情况进行调整。

运行代码得到的输出结果为方程的根值。

一、Newton法：#include<math.h>#include<stdio.h>double f(double x){return (3*x*x-exp(x));}double f1(double x){return (6*x-exp(x));}void main(){double x1=1,x;do{x=x1;x1=x-f(x)/f1(x);printf("x=%.9lf\n",x1);}while(fabs(x1-x)>0.000005);}说明：f 为原函数，f1为f的导函数，x1为初始值，通过"x=%.9lf“控制输入输出格式二、一般迭代法#include <stdio.h>#include <math.h>int main(){double x=1,x1;while(1){x1=pow(3*x+1,0.2);printf("x=%.6lf\n",x1);if(fabs(x1-x)<0.000005 )break;x=x1;}return 0;}说明：x1为初始值，x1=pow(3*x+1,0.2);为迭代格式，0.000005为允许误差，通过"x=%.6lf“控制输入输出格式三、Steffensen法：#include"stdio.h"#include"math.h"#define phi(x) pow(3*(x)+1,0.2);void main(){double x,x0,del,y,z;printf("x0="); scanf("%lf",&x0);printf("\ndel=:"); scanf("%lf",&del);while(1){y=phi(x0); z=phi(y);x=x0-(y-x0)*(y-x0)/(z-2*y+x0);printf("\n%.6lf",x);if(fabs(x-x0)<del) break;x0=x;}}说明：x0为初始值，pow(3*(x)+1,0.2);为φ（x）的格式，del为允许误差，通过"x=%.6lf“控制输入输出格式四、弦截法：#include<math.h>#include<stdio.h>double f(double x){ //计算f(x)的值return pow(2,x)+pow(3,x)-pow(4,x);}double point(double x1,double x2){//计算与x轴交点的x值printf("x=%.5f\n",(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1)));return (x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));}int main(){//输入两个数x1,x2double x1,x2,x;do{printf("输入两个数x1,x2:");scanf("%lf%lf",&x1,&x2);}while (f(x1)*f(x2)>= 0); // 当输入两个数大于0为真时，继续重新输入//关键循环步骤：do{x=point(x1,x2);//得到交点的值if(f(x)*f(x1)>0)x1=x;//新的x1elsex2=x;}while (fabs(f(x)) > 0.000005); }。

用c语言实现牛顿拉夫逊迭代法如何用C语言实现牛顿拉弗森迭代法[引言]牛顿拉弗森迭代法（Newton-Raphson method），又称为牛顿迭代法，是一种求解方程根的数值方法。

它是由17世纪英国物理学家伊萨克·牛顿和法国数学家约瑟夫·拉弗森独立提出的。

牛顿拉弗森迭代法通过不断逼近方程根，可以实现高精度的方程求解。

[原理]牛顿拉弗森迭代法的原理很简单，基于方程根附近的切线逼近曲线，通过不断迭代的方式逼近方程的根。

设方程的根为x0，代入方程得到曲线上的一个点(x0, f(x0))。

假设方程在x0附近可导，那么我们可以得到曲线在点(x0, f(x0))处的切线方程，即通过(x0, f(x0))和斜率f'(x0)得到切线的方程。

切线与x轴的交点即为新的近似根x1。

依此类推，我们可以采用迭代的方式不断逼近更精确的解。

具体迭代公式为：x_(n+1) = x_n - f(x_n) / f'(x_n)其中，x_n是第n次迭代得到的近似解，f(x_n)是方程在x_n处的函数值，f'(x_n)是方程在x_n处的导数值。

通过不断迭代，当两次迭代解的差值小于预设的误差范围时，我们可以认为得到了方程的近似根。

[步骤]接下来，让我们一步一步的使用C语言实现牛顿拉弗森迭代法。

步骤1：导入所需的头文件和函数声明c#include <stdio.h>#include <math.h>double f(double x); 定义方程的函数表达式double f_derivative(double x); 定义方程的导数函数表达式步骤2：实现方程和导数函数cdouble f(double x) {实现方程的函数表达式return x * x - 9;}double f_derivative(double x) {实现方程的导数函数表达式return 2 * x;}步骤3：实现牛顿拉弗森迭代法cdouble newton_raphson_method(double x0, double epsilon) { double x = x0;double difference = epsilon + 1; 初始化差值，使其大于预设的误差范围while (difference > epsilon) {double f_x = f(x);double f_derivative_x = f_derivative(x);x = x - f_x / f_derivative_x;difference = fabs(f(x));}return x;}步骤4：调用牛顿拉弗森迭代法求解方程cint main() {double x0 = 4; 初始值double epsilon = 0.0001; 误差范围double root = newton_raphson_method(x0, epsilon);printf("方程的近似根为：%lf\n", root);return 0;}[总结]通过以上步骤，我们可以使用C语言实现牛顿拉弗森迭代法来求解给定方程的根。

牛顿插值法的C 语言实现摘要：拉格朗日插值法具有明显的对称性，公式中的每一项与所有的插值节点有关。

因此，如果需要增加一个插值节点，则拉格朗日插值公式中的每一项都要改变，在有的应用中就显得不太方便。

因此，可以利用另外一种差值方法来弥补这种缺陷，就牛顿插值法。

本文通过对牛顿插值法的数学分析，主要给出其C 语言实现方法。

关键字：差商 差分 C 语言算法1差商及其牛顿插值公式1.1 差商及其主要性质定义 若已知函数()f x 在点(0,1,2,,)i x i n =⋅⋅⋅处的函数值()i f x 。

则称：00[]()f x f x =为函数()f x 在点0x 的0阶差商； 010101[][][,]f x f x f x x x x -=-为函数()f x 过点01,x x 的1阶差商；010201212[,][,][,,]f x x f x x f x x x x x -=-为函数()f x 过点012,,x x x 的2阶差商；以此类推，一般地称012101201211[,,,,][,,,,][,,,,,]k k k k k k k kf x x x x f x x x x f x x x x x x x -----⋅⋅⋅-⋅⋅⋅⋅⋅⋅=-为函数()f x 过点01,,x x ,⋅⋅⋅k x 的k 阶差商。

性质1 阶差商表示为函数值01(),(),,()k f x f x f x ⋅⋅⋅的线性组合。

即0120011()[,,,,]=()()()()kj k j j j j j j j k f x f x x x x x x x x x x x x =-+⋅⋅⋅-⋅⋅⋅--⋅⋅⋅-∑0()()kj kj ji i i jf x xx ==≠=-∑∏性质2 若函数()f x 在包含节点的区间[,]a b 上存在k 阶导数，则k 阶差商与导数的关系为()01()[,,,],[,]!k k f f x x x a b k ξξ⋅⋅⋅=∈ 1.2 牛顿插值公式通过1n +个互异点上的次数不超过n 的插值多项式()n P x 可以表示如下形式：0010101201()[][,]()[,,]()()n P x f x f x x x x f x x x x x x x =+-+--+⋅⋅⋅0101[,,,]()()n n f x x x x x x x -+⋅⋅⋅-⋅⋅⋅-这种形式的插值多项式称为牛顿插值多项式，一般记为0010101201()[][,]()[,,]()()n N x f x f x x x x f x x x x x x x =+-+--+⋅⋅⋅0101[,,,]()()n n f x x x x x x x -+⋅⋅⋅-⋅⋅⋅-由牛顿插值多项式可以看出，当增加一个插值点时，当前已有的各项不变，只需要在后面增加一项即可。

用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();}。

简介：本程序用牛顿插值法对函数表，X值选取为1.5测试本程序。

源程序：#include <iostream.h>#include <math.h>void main(){int I;char L;do{double M[100][100];double x[100],y[100];double X=1,xx=0,w=1,N=0,P,R=1;int n;cout<<"请输入所求均差阶数：";cin>>n;for(int i=0;i<=n;i++){cout<<"请输入x"<<i<<"的值："<<endl; cin>>x[i];cout<<"请输入y"<<i<<"的值："<<endl; cin>>y[i];M[i][0]=x[i];M[i][1]=y[i];}for( int j=2;j<=n+1;j++){for( i=1;i<=n;i++){M[i][j]=(M[i][j-1]-M[i-1][j-1])/(M[i][0]-M[i-j+1][0]);}}for(i=1;i<=n;i++){cout<<"其"<<i<<"阶均差为："<<M[i][i+1]<<endl;}cout<<"请输入x的值：x=";cin>>xx;for(i=0;i<n;i++){X*=xx-x[i];N+=M[i+1][i+2]*X;P=M[0][1]+N;}cout<<"其函数值：y="<<P<<endl;cout<<endl<<"如还想算其它插值请按'y'否则按'n'"<<endl; cin>>L;}while(L=='y');}测试结果：。

⽜顿插值法及其C++实现⽜顿插值法⼀、背景引⼊相信朋友们，开了拉格朗⽇插值法后会被数学家的思维所折服，但是我想说有了拉格朗⽇插值法还不够，因为我们每次增加⼀个点都得重算所有插值基底函数，这样会增加计算量，下⾯我们引⼊⽜顿插值法，这种插值法，添加⼀个插值结点我们只要做很⼩的变动便可以得到新的插值多项式。

⼆、理论推导-均差的定义：（⼀阶均差)⼆阶均差为⼀阶均差再求均差。

（显然是递推的)⼀般地，函数f 的k阶均差定义为：由均差的性质可以推导出：k+1阶均差:(具体性质看:《数值分析：第5版》 page:30)由均差的递推性，我们可以⽤以下表来求:求表的公式：table[i][j] = (table[i - 1][j] - table[i - 1][j - 1]) / (x[j] - x[j - i]);其中P（x) 为插值多项式，⽽R（x) 为插值余项。

所以p（x）：（由于图⽚问题此处P（x) 同N(x))三、代码实现由以上推导可知，求⽜顿插值多项式⼦主要就是求均差。

均差可由上表递推求得：求表的公式：table[i][j] = (table[i - 1][j] - table[i - 1][j - 1]) / (x[j] - x[j - i]);#include <iostream>using namespace std;#include <vector>inline double newton_solution(double x[], double y[], int n, double num, int newton_time) {vector<vector<double> > table(n + 1);for (int i = 0; i <= n; i++) {table[i].resize(n + 1);}for (int i = 0; i <= n; i++) table[0][i] = y[i];for (int i = 1; i <= n; i++) {for (int j = i; j <= n; j++) {table[i][j] = (table[i - 1][j] - table[i - 1][j - 1]) / (x[j] - x[j - i]); }}double res = 0.0;for (int i = 0; i <= newton_time; i++) {double temp = table[i][i];for (int j = 0; j < i; j++) {temp *= num - x[j];}res += temp;}return res;}int main(int argc, char const *argv[]){int n = 0;cout << "插值节点个数-1:";cin >> n;double x[n + 1], y[n + 1];cout << "\n请输⼊x[i]:";for (int i = 0; i <= n; i++) {cin >> x[i];}cout << "\n请输⼊y[i]:";for (int i = 0; i <= n; i++) {cin >> y[i];}double num = 0;cout << "\n请输⼊要求的点的x:";cin >> num;cout << "\n请输⼊所求的插值多项式次数:";double newton_time = 0;cin >> newton_time;cout << newton_solution(x, y, n, num, newton_time) << endl; return0;。

计算方法与实习上机实验（三）实验目的：编写一个用牛顿前插公式计算函数值的程序，要求先输出差分表，再计算x点的函数值，并应用于下面的问题：算法:(1).输入n=4,xi,yi(i=0,1,2,3,4).(2).计算各阶差分f00,f10,f20,f30,f40。

(3)计算f00+f10*t+f20*t*(t-1)/2+f30*t*(t-1)*(t-2)/6+f40*t*(t-1)*(t-2)*(t-3)/2 实验程序：#include<iostream>using namespace std;double NewTon(float f00,float f10,float f20,float f30,float f40,float t){return f00+f10*t+f20*t*(t-1)/2+f30*t*(t-1)*(t-2)/6+f40*t*(t-1)*(t-2)*(t-3)/24;}void main(){ double f1[4],f2[3],f3[2],f4[1];double f0[5]={1.30103,1.32222,1.34242,1.36173,1.38021};double x[5]={20,21,22,23,24};cout<<" "<<"差分表"<<endl;for(int i=0;i<=4;i++){cout<<" x"<<i<<"="<<x[i];}cout<<endl<<" ";for(i=0;i<=4;i++){cout<<" "<<f0[i];}cout<<endl<<"1阶差分"<<" ";for(i=0;i<=3;i++){ f1[i]=f0[i+1]-f0[i];cout<<f1[i]<<" ";}cout<<endl<<"2阶差分"<<" ";for(i=0;i<=2;i++){ f2[i]=f1[i+1]-f1[i];cout<<f2[i]<<" ";}cout<<endl<<"3阶差分"<<" ";for(i=0;i<=1;i++){ f3[i]=f2[i+1]-f2[i];cout<<f3[i]<<" ";}cout<<endl<<"4阶差分"<<" ";for(i=0;i<1;i++){ f4[i]=f3[i+1]-f3[i];cout<<f4[i]<<endl;}double t=21.4-x[0];cout<<"N(21.4)="<<NewTon(f0[0],f1[0],f2[0],f3[0],f4[0],t)<<endl; }实验结果：。

计算机数值计算方法及程序设计P63函数值为//拉格朗日插值P63#include <stdio.h>float czl(int n,float x1,float *px,float *py);int main(){float x1,y1;int n;float *p1,*p2;float x[10]={1,2,3,4,5,6,7};float y[10]={1,1.414214,1.732051,2,2.236068,2.449490,2.645751};printf("Input numbers:x1 n=\n");scanf("%f%d",&x1,&n);p1=x;p2=y;y1=czl(n,x1,p1,p2);printf("y1=%f\n",y1);getch();return 0;}float czl(int n,float x1,float *px,float *py){int i,j;float x[10],y[10],t,y1;y1=0.0;for(i=0;i<n;i++,px++,py++){x[i]=*px;y[i]=*py;}for(i=0;i<n;i++){t=1.0;for(j=0;j<n;j++)if(i!=j)t=t*(x1-x[j])/(x[i]-x[j]);y1=y1+t*y[i];}return(y1);}//输入为//2.5 4//输出为//y=1.582274//线性插值P58#include <stdio.h>float cz(float x0,float x1,float y0,float y1,float x); int main(void){float x0,x1,y0,y1,x,y;printf("Input numbers:x0,x1,y0,y1,x=?\n");scanf("%f%f%f%f%f",&x0,&x1,&y0,&y1,&x);y=cz(x0,x1,y0,y1,x);printf("y=%f\n",y);}float cz(float x0,float x1,float y0,float y1,float x) {float l0,l1,y;l0=(x-x1)/(x0-x1);l1=(x-x0)/(x1-x0);y=l0*y0+l1*y1;return (y);}/*输入：1 5 1 2.6.68 2.5输出y=1.463526 *////抛物插值P60#include<stdio.h>float cz(float x0,float x1,float x2,float y0,float y1,float y2,float x); float cz(float x0,float x1,float x2,float y0,float y1,float y2,float x) {float l0,l1,l2,y;l0=(x-x1)*(x-x2)/((x0-x1)*(x0-x2));l1=(x-x0)*(x-x2)/((x1-x0)*(x1-x2));l2=(x-x0)*(x-x1)/((x2-x0)*(x2-x1));y=l0*y0+l1*y1+l2*y2;return(y);}int main(void){float x0,x1,x2,y0,y1,y2,x,y;printf("Input numbers:x0 x1 x2 y0 y1 y2 x=?\n");freopen("in.txt","r",stdin);freopen("out.txt","w",stdout);scanf("%f%f%f%f%f%f%f",&x0,&x1,&x2,&y0,&y1,&y2,&x);y=cz(x0,x1,x2,y0,y1,y2,x);printf("y=%f\n",y);getch();getch();return 0;}/*输入：1 3 5 1 1.732051 2.236068 2.5输出y=1.570416 *///牛顿插值P83#include<stdio.h>#include<math.h>#define N 6float sub(float a[],float b[],float x,float e); main(){float u[N]={100,121,122,169,196,225}; float v[N]={10,11,12,13,14,15};float x,y,e,*p1,*p2;printf("Input number x e= ");scanf("%f%e",&x,&e);p1=u;p2=v;y=sub(p1,p2,x,e);printf("y=%f\n",y);getch();getch();}float sub(float *pp1,float *pp2,float x,float e) {float a[N],b[N],t[N],y,y1,c;int i,k;for(i=0;i<N;i++,pp1++) {a[i]=*pp1;printf("%12.6",a[i]);}printf("\n");for(i=0;i<N;i++,*pp2++) {b[i]=*pp2;printf("%12.6f",b[i]);}printf("\n");y1=b[0];y=0;c=1.0;for(k=1;k<N;k++) {for(i=k;i<N;i++){t[i]=(b[i]-b[i-1])/(a[i]-a[i-k]);}c=c*(x-a[k-1]);y1=y1+c*t[k];if(fabs(y-y1)<e) y=y1;for(i=k;i<N;i++){b[i]=t[i];}}return(y);}。