东南大学计算方法与实习上机实验三

计算方法B上机实习报告计算方法B 上机实习报告1.对以下和式计算：0142118184858616n n S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求： (1)若只需保留11个有效数字，该如何进行计算；(2)若要保留30个有效数字，则又将如何进行计算；问题分析：在该题中S 的每一项14211()1681848586n ns n n n n =---++++存在两个相近的数相减的问题，因此为了避免有效数字损失，最好是改变运算顺序，分别将正数和负数相加，然后再将其和相加。

另外，sn 中有多个负数相加，可以按照绝对值递增的顺序求和，以减少舍入误差的影响。

同时，为了避免大数吃小数的问题，本题先计算出保留目标有效数字所需 要的迭代次数，然后采用倒序相加的方法实现。

程序实现：clear;clc;m=input('请输入要保留的有效数字位数：'); s1=0; s2=0; k=0; s=1;%%%%判断多需要的迭代次数 while s>=0.5*10^-(m-1)s=4/(16^k*(8*k+1))-(2/(16^k*(8*k+4))+1/(16^k*(8*k+5))+1/(16^k*(8*k+6))); k=k+1; end%%%%正负数分别按照绝对值递增的顺序倒序相加 for n=(k-1):-1:0 a1=4/(16^n*(8*n+1)); a2=2/(16^n*(8*n+4)); a3=1/(16^n*(8*n+5)); a4=1/(16^n*(8*n+6)); s1=a1+s1; s2=a4+a3+a2+s2; end S=s1-s2; S=vpa(S,m)运算结果：总结心得：在计算求和问题中，应特别注意相近数相减的问题，这样会造成有效数字灾难性的损失。

另外在两个数量级相差较大的数字相加减时，较小数的有效数字会被丧失，因此要按照从小到大的顺序相加。

在上题计算中分别对正负相采用倒序相加，这样就有效的避免了“大数吃小数”的问题。

计算方法A上机报告学院（系）：电气工程学院学生姓名：陶然学号：授课老师：完成日期：2019年12月03日西安交通大学Xi'an Jiaotong University目录1 QR分解法求解线性方程组 (2)1.1 算法原理 (2)1.1.1 基于吉文斯变换的QR分解 (2)1.1.2 基于豪斯霍尔德变换的QR分解 (3)1.2 程序流程图 (4)1.2.1 基于吉文斯变换的QR分解流程图 (4)1.2.2 基于豪斯霍尔德变换的QR分解流程图 (5)1.3 程序使用说明 (5)1.3.1 基于吉文斯变换的QR分解程序说明 (5)1.3.2 基于豪斯霍尔德变换的QR分解程序说明 (7)1.4 算例计算结果 (8)2 共轭梯度法求解线性方程组 (10)2.1 算法原理 (10)2.2 程序流程图 (10)2.3 程序使用说明 (11)2.4 算例计算结果 (12)3 三次样条插值 (14)3.1 算法原理 (14)3.2 程序流程图 (16)3.3 程序使用说明 (17)3.4 算例计算结果 (19)4 龙贝格积分 (21)4.1 算法原理 (21)4.2 程序流程图 (22)4.3 程序使用说明 (23)4.4 算例计算结果 (24)结论 (26)1 QR 分解法求解线性方程组1.1 算法原理矩阵的QR 分解是指，可以将矩阵A 分解为一个正交矩阵Q 和一个上三角矩阵R 的乘积，实际中，QR 分解经常被用来解决线性最小二乘问题，分解情况如图1.1所示。

=⨯图1.1 QR 分解示意图本次上机学习主要进行了两个最基本的正交变换—吉文斯（Givens ）变换和豪斯霍尔德（Householder ）变换，并由此导出矩阵的QR 分解以及求解线性方程组的的方法。

1.1.1 基于吉文斯变换的QR 分解吉文斯矩阵也称初等旋转阵，如式（1.1）所示，它把n 阶单位矩阵I 的第,i j 行的对角元改为c ，将第i 行第j 列的元素改为s ，第j 行第i 列的元素改为s −后形成的矩阵。

计算方法A上机实验报告姓名：苏福班级：硕4020 学号：3114161019一、上机练习目的1）复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。

2）利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。

二、上机练习任务1）利用计算机语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。

2）掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。

3）写出上机练习报告。

三、上机题目1. 共轭梯度法求解线性方程组。

（第三章）2. 三次样条插值（第四章）3. 龙贝格积分（第六章）4. 四阶龙格－库塔法求解常微分方程的初值问题四、上机报告题目1：共轭梯度法求解线性方程组1．算法原理共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小值的问题。

从任意给定的初始点出发，沿一组关于矩阵A共轭的方向进行线性搜索，在无舍入误差的假定下，最多迭代n 次（其中n 为矩阵A 的阶数），就可求得二次函数的极小值，也就求得了线性方程组Ax b =的解。

定理：设A 是n 阶对称正定矩阵，则x *是方程组Ax b =的解得充分必要条件是x *是二次函数1()2TT f x x Ax b x =-的极小点，即 ()()min nx R Ax b f x f x **∈=⇔=共轭梯度法的计算公式：(0)(0)(0)()()()()(1)()()(1)(1)(1)()()()(1)(1)()k T k k k T k k k k k k k k T k k k T k k k k k d r b Ax r d d Ad xx d r b Ax r Ad d Ad d r d ααββ++++++⎧==-⎪⎪=⎪⎪=+⎪⎨=-⎪⎪⎪=-⎪⎪=+⎩2. 程序框图（1）编写共轭梯度法求解对称正定矩阵的线性方程组见附录（myge.m）：function x=myge(A,b)输入对称正定矩阵及对应的列向量，初始向量设为0，精度取为810 。

1.拉格朗日插值多项式，用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>float lagrange(float *x,float *y,float xx,int n) /*拉格朗日插值算法*/{ int i,j;float *a,yy=0.0; /*a作为临时变量，记录拉格朗日插值多项式*/a=(float *)malloc(n*sizeof(float));for(i=0;i<=n-1;i++){ a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}free(a);return yy;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error!The value of n must in (0,20)."); getch();return 1;} if(n<=0) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");printf("Input xx:");scanf("%f",&xx);yy=lagrange(x,y,xx,n);printf("x=%f,y=%f\n",xx,yy);getch();}2.牛顿插值多项式，用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>void difference(float *x,float *y,int n){ float *f;int k,i;f=(float *)malloc(n*sizeof(float));for(k=1;k<=n;k++){ f[0]=y[k];for(i=0;i<k;i++)f[i+1]=(f[i]-y[i])/(x[k]-x[i]);y[k]=f[k];}return;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} if(n<=0) {printf("Error! The value of n must in (0,20).");getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");difference(x,(float *)y,n);printf("Input xx:");scanf("%f",&xx);yy=y[20];for(i=n-1;i>=0;i--) yy=yy*(xx-x[i])+y[i];printf("NewtonInter(%f)=%f",xx,yy);getch();}3.高斯列主元消去法,求解其次线性方程组第一种#include<stdio.h>#include <math.h>#define N 20int main(){ int n,i,j,k;int mi,tmp,mx;float a[N][N],b[N],x[N];printf("\nInput n:");scanf("%d",&n);if(n>N){ printf("The input n should in(0,N)!\n");getch();return 1;}if(n<=0){ printf("The input n should in(0,N)!\n");getch();return 1;}printf("Now input a(i,j),i,j=0...%d:\n",n-1); for(i=0;i<n;i++){ for(j=0;j<n;j++)scanf("%f",&a[i][j]);}printf("Now input b(i),i,j=0...%d:\n",n-1);for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n-2;i++){ for(j=i+1,mi=i,mx=fabs(a[i][j]);j<n-1;j++) if(fabs(a[j][i])>mx){ mi=j;mx=fabs(a[j][i]);}if(i<mi){ tmp=b[i];b[i]=b[mi];b[mi]=tmp;for(j=i;j<n;j++){ tmp=a[i][j];a[i][j]=a[mi][j];a[mi][j]=tmp;}}for(j=i+1;j<n;j++){ tmp=-a[j][i]/a[i][i];b[j]+=b[i]*tmp;for(k=i;k<n;k++)a[j][k]+=a[i][k]*tmp;}}x[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i>=0;i--){ x[i]=b[i];for(j=i+1;j<n;j++)x[i]-=a[i][j]*x[j];x[i]/=a[i][i];}for(i=0;i<n;i++)printf("Answer:\n x[%d]=%f\n",i,x[i]); getch();return 0;}第二种#include<math.h>#include<stdio.h>#define NUMBER 20#define Esc 0x1b#define Enter 0x0dfloat A[NUMBER][NUMBER+1] ,ark;int flag,n;exchange(int r,int k);float max(int k);message();main(){float x[NUMBER];int r,k,i,j;char celect;clrscr();printf("\n\nUse Gauss.");printf("\n\n1.Jie please press Enter."); printf("\n\n2.Exit press Esc.");celect=getch();if(celect==Esc)exit(0);printf("\n\n input n=");scanf("%d",&n);printf(" \n\nInput matrix A and B:"); for(i=1;i<=n;i++){printf("\n\nInput a%d1--a%d%d and b%d:",i,i,n,i);for(j=1;j<=n+1;j++) scanf("%f",&A[i][j]);}for(k=1;k<=n-1;k++){ark=max(k);if(ark==0){printf("\n\nIt's wrong!");message();}else if(flag!=k)exchange(flag,k);for(i=k+1;i<=n;i++)for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k];}x[n]=A[n][n+1]/A[n][n];for( k=n-1;k>=1;k--){float me=0;for(j=k+1;j<=n;j++){me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}for(i=1;i<=n;i++){printf(" \n\nx%d=%f",i,x[i]);}message();}exchange(int r,int k){int i;for(i=1;i<=n+1;i++)A[0][i]=A[r][i];for(i=1;i<=n+1;i++)A[r][i]=A[k][i];for(i=1;i<=n+1;i++)A[k][i]=A[0][i];}float max(int k){int i;float temp=0;for(i=k;i<=n;i++)if(fabs(A[i][k])>temp){temp=fabs(A[i][k]);flag=i;}return temp;}message(){printf("\n\n Go on Enter ,Exit press Esc!");switch(getch()){case Enter: main();case Esc: exit(0);default:{printf("\n\nInput error!");message();} }}4.龙贝格求积公式，求解定积分#include<stdio.h>#include<math.h>#define f(x) (sin(x)/x)#define N 20#define MAX 20#define a 2#define b 4#define e 0.00001float LBG(float p,float q,int n){ int i;float sum=0,h=(q-p)/n;for (i=1;i<n;i++)sum+=f(p+i*h);sum+=(f(p)+f(q))/2;return(h*sum);}void main(){ int i;int n=N,m=0;float T[MAX+1][2];T[0][1]=LBG(a,b,n);n*=2;for(m=1;m<MAX;m++){ for(i=0;i<m;i++)T[i][0]=T[i][1];T[0][1]=LBG(a,b,n);n*=2;for(i=1;i<=m;i++)T[i][1]=T[i-1][1]+(T[i-1][1]-T[i-1][0])/(pow(2,2*m)-1);if((T[m-1][1]<T[m][1]+e)&&(T[m-1][1]>T[m][1]-e)){ printf("Answer=%f\n",T[m][1]); getch();return ;}}}5.牛顿迭代公式，求方程的近似解#include<stdio.h>#include<math.h>#include<conio.h>#define N 100#define PS 1e-5#define TA 1e-5float Newton(float (*f)(float),float(*f1)(float),float x0 ) { float x1,d=0;int k=0;do{ x1= x0-f(x0)/f1(x0);if((k++>N)||(fabs(f1(x1))<PS)){ printf("\nFailed!");getch();exit();}d=(fabs(x1)<1?x1-x0:(x1-x0)/x1);x0=x1;printf("x(%d)=%f\n",k,x0);}while((fabs(d))>PS&&fabs(f(x1))>TA) ;return x1;}float f(float x){ return x*x*x+x*x-3*x-3; }float f1(float x){ return 3.0*x*x+2*x-3; }void main(){ float f(float);float f1(float);float x0,y0;printf("Input x0: ");scanf("%f",&x0);printf("x(0)=%f\n",x0);y0=Newton(f,f1,x0);printf("\nThe root is x=%f\n",y0); getch();}。

东南大学自动化学院实验报告课程名称：信息通信网络概论第3次实验实验名称：实验三基于客户/服务器模式的网络通信编程实现院（系）：自动化专业：自动化姓名：学号：实验室：金智楼实验组别：同组人员：实验时间： 2016 年 12 月 13 日评定成绩：审阅教师：目录一．实验目的和要求 (3)二．实验原理 (3)三. 实验方案与实验步骤 (4)四．实验设备与器材配置 (5)五．实验记录 (5)六．实验总结 (10)七．思考题或讨论题 (11)附录：部分代码一．实验目的和要求1.进一步了解网络编程的过程；2.掌握Windows环境下基于WinSock的编程方法和通信实现；3.熟悉客户/服务器模式的网络通信编程实现，编写一个聊天工具，即以客户端和服务器端的模式进行互发消息。

二．实验原理一个在建立分布式应用时最常用的范例便是客户机/服务器模型。

在这种方案中客户应用程序向服务器程序请求服务。

这种方式隐含了在建立客户机/服务器间通讯时的非对称性。

客户机/服务器模型工作时要求有一套为客户机和服务器所共识的惯例来保证服务能够被提供（或被接受）。

这一套惯例包含了一套协议。

它必须在通讯的两头都被实现。

根据不同的实际情况，协议可能是对称的或是非对称的。

在对称的协议中，每一方都有可能扮演主从角色；在非对称协议中，一方被不可改变地认为是主机，而另一方则是从机。

一个对称协议的例子是Internet中用于终端仿真的TELNET。

而非对称协议的例子是Internet中的FTP。

无论具体的协议是对称的或是非对称的，当服务被提供时必然存在“客户进程”和“服务进程”。

一个服务程序通常在一个众所周知的地址监听对服务的请求，也就是说，服务进程一直处于休眠状态，直到一个客户对这个服务的地址提出了连接请求。

在这个时刻，服务程序被“惊醒”并且为客户提供服务－对客户的请求作出适当的反应。

这一请求/相应的过程可以简单的用图2-1表示。

虽然基于连接的服务是设计客户机/服务器应用程序时的标准，但有些服务也是可以通过数据报套接口提供的。

《数值计算方法》上机实验报告华北电力大学实验名称数值il•算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一.各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程*对于非线性方程，若已知根的一个近似值，将在处展开成一阶xxfx ()0, fx ()xkk泰勒公式"f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2!忽略高次项，有，fxfxfxxx 0 ()()()，，, kkk右端是直线方程，用这个直线方程来近似非线性方程。

将非线性方程的**根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkkfx 0 fx 0 0,解出fX 0 *k XX,, k' fx 0 k水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ikfx ()k 八XX, Ikk* fx()k这就是牛顿迭代公式。

,2,计算机程序框图:，见,,3,输入变量、输出变量说明:X输入变量:迭代初值，迭代精度，迭代最大次数，\0输出变量:当前迭代次数，当前迭代值xkl,4,具体算例及求解结果:2/16华北电力大学实验报吿开始读入l>k/fx()0?,0fx 0 Oxx，，01* fx ()0XX,,，?10kk, ,1，kN, ?xx, 10输出迭代输出X输出奇异标志1失败标志,3,输入变量、输出变量说明: 结束例：导出计算的牛顿迭代公式，并il •算。

（课本P39例2-16） 115cc （0）, 求解结果:10. 75000010.72383710. 72380510. 7238052、列主元素消去法求解线性方程组，1,算法原理:高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角3/16华北电力大学实验报告方程组求解。

实习题11. 用两种不同的顺序计算644834.11000012≈∑=-n n，试分析其误差的变化解：从n=1开始累加，n 逐步增大，直到n=10000；从n=10000开始累加，n 逐步减小，直至1。

算法1的C 语言程序如下： #include<stdio.h> #include<math.h> void main() { float n=0.0; int i; for(i=1;i<=10000;i++) { n=n+1.0/(i*i); } printf("%-100f",n); printf("\n"); float m=0.0; int j; for(j=10000;j>=1;j--) { m=m+1.0/(j*j); } printf("%-7f",m); printf("\n"); }运行后结果如下：结论： 4.设∑=-=Nj N j S 2211，已知其精确值为)11123(21+--N N 。

1）编制按从大到小的顺序计算N S 的程序； 2）编制按从小到大的顺序计算N S 的程序；3）按2种顺序分别计算30000100001000,,S S S ，并指出有效位数。

解：1）从大到小的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=N;i>1;i--){n=n+1.0/(i*i-1);N=N-1;}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：2）从小到大的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=2;i<=N;i++){n=n+1.0/(i*i-1);}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：结论：通过比较可知：N 的值较小时两种算法的运算结果相差不大，但随着N 的逐渐增大，两种算法的运行结果相差越来越大。

计算方法与实习上机实验报告一、引言本文旨在介绍和展示我们在“计算方法”课程中的实习上机实验环节所完成的一些关键任务和所取得的成果。

该实验课程的目标是让我们更深入地理解和应用各种计算方法，并在实际操作中提高我们的编程和问题解决能力。

二、实验内容与目标实验的主要内容是利用各种计算方法解决实际数学问题。

我们被要求使用编程语言（如Python或Java）来实现和解决这些问题。

这些问题包括使用牛顿法求解平方根，使用蒙特卡洛方法计算圆周率，以及使用最优化方法求解函数的最小值等。

实验的目标不仅是让我们掌握计算方法的基本理论，更是要让我们能够在实际操作中运用这些方法。

我们需要在实习过程中，通过与同伴们合作，共同解决问题，提高我们的团队合作能力和问题解决能力。

三、实验过程与问题解决策略在实验过程中，我们遇到了许多问题，如编程错误、理解困难和时间压力等。

我们通过相互讨论、查阅资料和寻求教师帮助等方式，成功地解决了这些问题。

例如，在实现牛顿法求解平方根时，我们一开始对导数的计算和理解出现了一些错误。

但我们通过查阅相关资料和讨论，最终理解了导数的正确计算方法，并成功地实现了牛顿法。

四、实验结果与结论通过这次实习上机实验，我们不仅深入理解了计算方法的基本理论，还在实际操作中提高了我们的编程和问题解决能力。

我们的成果包括编写出了能有效求解平方根、计算圆周率和求解函数最小值的程序。

这次实习上机实验非常成功。

我们的团队不仅在理论学习和实践操作上取得了显著的进步，还在团队合作和问题解决方面积累了宝贵的经验。

这次实验使我们对计算方法有了更深的理解和认识，也提高了我们的编程技能和解决问题的能力。

五、反思与展望回顾这次实验，我们意识到在实验过程中，我们需要更好地管理我们的时间和压力。

在解决问题时，我们需要更有效地利用我们的知识和资源。

在未来，我们希望能够更加熟练地运用计算方法，并能够更有效地解决问题。

我们也希望能够将所学的计算方法应用到更广泛的领域中，如数据分析、科学研究和工业生产等。

计算方法与实习实验报告目录实习题1-1 (2)实习题2-1（1） (3)实习题3-2 (4)实习题3-5（1） (6)实习题4-2 (9)实习题5-1 (10)实习题6 (13)实习题1-1用2种不同的顺序计算644834.11000012≈∑=-n n，分析其误差的变化。

算法1:从n=1开始累加，n 逐步增大，直至n=10000 算法2:从n=10000开始累加，n 逐步减小，直至n=1 算法1的Matlab 程序如下： m=0;for n=1:10000 yl=1/(n.*n); m=m+yl; endformat long m算法1的输出结果m=1.644834071848065算法2的Matlab 程序如下： m=0;for n=10000:-1:1 yl=1/(n.*n); m=m+yl; endformat long m算法2的输出结果m=1.644834071848060结论分析：由运行结果可以发现，当从n=1开始计算时，算出结果为1.644834071848065，而当从n=10000开始计算时，算出结果为1.644834071848060.与题中所给结果相比，可知从n=10000开始计算时，所得近似结果更接近真实值。

原因是从n=1开始算时，由于被累加的数是越来越小的，而且开始的几个数比后面的数大得多，而浮点数保留位数有限，因此在累加过程中出现了大数吃小数的现象，导致误差较大。

而从n=10000开始累加就很好地避免以上问题，从而能够更接近精确值。

实习题2-1（1）用牛顿法求方程02=-x e x 的根算法：给定初始值0x ,ε为根的容许误差，η为)(x f 的容许误差，N 为迭代次数的容许值。

①如果0)(0'=x f 或迭代次数大于N ，则算法失败，结束；否则执行②。

②计算)()(0'001x f x f x x -=。

③若01x x -<ε或)(1x f <η，则输出1x ,程序结束；否则执行④。

数值分析上机实验报告目录1.chapter1舍入误差及有效数 (1)2.chapter2Newton迭代法 (3)3.chapter3线性代数方程组数值解法-列主元Gauss消去法 (7)4.chapter3线性代数方程组数值解法-逐次超松弛迭代法 (8)5.chapter4多项式插值与函数最佳逼近 (10)1.chapter1舍入误差及有效数1.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。

（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。

（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。

（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。

（编制程序时用单精度）（4）通过本次上机题，你明白了什么？ 1.2编写相应的matlab 程序 clear;N=input('please input N:'); AValue=((3/2-1/N-1/(N+1))/2); sn1=single(0); sn2=single(0); for i=2:Nsn1=sn1+1/(i*i-1); %从大到小相加的通用程序% endep1=abs(sn1-AValue); for j=N:-1:2sn2=sn2+1/(j*j-1); %从小到大相加的通用程序% endep2=abs(sn2-AValue);fprintf('精确值为:%f\n',AValue);fprintf('从大到小的顺序累加得sn=%f\n',sn1); fprintf('从大到小相加的误差ep1=%f\n',ep1); fprintf('从小到大的顺序累加得sn=%f\n',sn2); fprintf('从小到大相加的误差ep2=%f\n',ep2); disp('================================='); 1.3matlab 运行程序结果 >> chaper1please input N:100 精确值为:0.740050从大到小的顺序累加得sn=0.740049 从大到小相加的误差ep1=0.000001 从小到大的顺序累加得sn=0.740050 从小到大相加的误差ep2=0.000000 >> chaper1please input N:10000 精确值为:0.749900从大到小的顺序累加得sn=0.749852 从大到小相加的误差ep1=0.000048 从小到大的顺序累加得sn=0.749900 从小到大相加的误差ep2=0.000000please input N:1000000精确值为:0.749999从大到小的顺序累加得sn=0.749852 从大到小相加的误差ep1=0.000147 从小到大的顺序累加得sn=0.749999 从小到大相加的误差ep2=0.0000001.4结果分析以及感悟按照从大到小顺序相加的有效位数为：5,4,3。

. 《计算方法》上机实验报告班级：XXXXXX小组成员：XXXXXXXXXXXXXXXXXXXXXXXXXXXX任课教师：XXX二〇一八年五月二十五日前言通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。

以下为本次上机实验报告，按照实验内容共分为六部分。

实验一：一、实验名称及题目： Newton 迭代法例2.7(P38):应用Newton 迭代法求在附近的数值解，并使其满足.二、解题思路：设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标)(')(0001x f x f x x -=,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标)(')(1112x f x f x x -=称2x 为'x 的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把)(')(1n n n n x f x f x x -=+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。

三、Matlab 程序代码：function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0);y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1;while abs(x1-x0)>=tol x0=x1;y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; endx=double(x1) K四、运行结果：实验二：一、实验名称及题目：Jacobi 迭代法例3.7（P74）：试利用Jacobi 迭代公式求解方程组要求数值解为方程组的精确解. 二、解题思路：首先将方程组中的系数矩阵A 分解成三部分，即：U D L A ++=，D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。

东南大学自动化学院实验报告课程名称：计算机控制技术第三次实验实验名称：离散化方法研究院（系）：自动化专业：自动化姓名：学号：同组人员：实验时间：2017 年 4 月12 日评定成绩：审阅教师：目录一．实验目的 (3)二．实验设备 (3)三．实验原理 (3)四．实验步骤 (7)五．实验结果 (8)一、实验目的1．学习并掌握数字控制器的设计方法（按模拟系统设计方法与按离散设计方法）；2．熟悉将模拟控制器D(S)离散为数字控制器的原理与方法（按模拟系统设计方法）；3．通过数模混合实验，对D(S)的多种离散化方法作比较研究，并对D(S)离散化前后闭环系统的性能进行比较，以加深对计算机控制系统的理解。

二、实验设备1．THBDC-1型控制理论·计算机控制技术实验平台2．PCI-1711数据采集卡一块3．PC机1台(安装软件“VC++”及“THJK_Server”)三、实验原理由于计算机的发展，计算机及其相应的信号变换装置（A/D和D/A）取代了常规的模拟控制。

在对原有的连续控制系统进行改造时，最方便的办法是将原来的模拟控制器离散化。

在介绍设计方法之前，首先应该分析计算机控制系统的特点。

图3-1为计算机控制系统的原理框图。

图3-1 计算机控制系统原理框图由图3-1可见，从虚线I向左看，数字计算机的作用是一个数字控制器，其输入量和输出量都是离散的数字量，所以，这一系统具有离散系统的特性，分析的工具是z变换。

由虚线II向右看，被控对象的输入和输出都是模拟量，所以该系统是连续变化的模拟系统，可以用拉氏变换进行分析。

通过上面的分析可知，计算机控制系统实际上是一个混合系统，既可以在一定条件下近似地把它看成模拟系统，用连续变化的模拟系统的分析工具进行动态分析和设计，再将设计结果转变成数字计算机的控制算法。

也可以把计算机控制系统经过适当变换，变成纯粹的离散系统，用z变化等工具进行分析设计，直接设计出控制算法。

按模拟系统设计方法进行设计的基本思想是，当采样系统的采样频率足够高时，采样系统的特性接近于连续变化的模拟系统，此时忽略采样开关和保持器，将整个系统看成是连续变化的模拟系统，用s 域的方法设计校正装置D(s)，再用s 域到z 域的离散化方法求得离散传递函数D(z)。

东南大学计算方法与实习实验报告计算方法与实习实验报告学院:学号:姓名：完成日期:实习题一4、设2211Nn j S j ==-∑，已知其精确值为。

1）编制按从大到小的顺序计算S n 的程序； 2）编制按从小到大的顺序计算S n 的程序；3）按两种顺序分别计算S 1000，S 10000，S 30000，并指出有效位数。

● 实验代码 C 语言程序如下：#include #include using namespace std; int main(){ float Sn=0; int N; cin>>N; for(float j=2;j<=N;j++){ Sn=1/(j*j-1)+Sn; } cout<<"从小到大计算的结果为"<<sn<for(j=N;j>=2;j--){ Sn=1/(j*j-1)+Sn;}cout<<"从大到小计算的结果为"<<sn<<=""> ● 运行窗口实习题二1、用牛顿法求下列方程的根：1） 20xx e -=实验代码C 语言程序代码如下：#include #include #define N 100 #define eps 1e-6 #define eta 1e-8using namespace std;float Newton(float f(float),float fl(float),float x0){ float x1,d; int k=0;do{x1=x0-f(x0)/fl(x0);if(k++>N||fabs(fl(x1))<eps){cout<<"发散"<<endl;break;}< p="">d=fabs(x1)<1?x1-x0:(x1-x0)/x1;x0=x1;cout<<"x="<<x0<<endl;< p="">}while(fabs(d)>eps&&fabs(f(x1))>eta);return x1;}float f(float x){return x+log10(x)-2;}float fl(float x){return 1+1/x;}void main(){float x0,y0;cin>>x0;y0=Newton(f,fl,x0);cout<<"方程的根为"<<y0<<endl;< p="">}运行窗口实习题三1、用列主元消去法解方程组：1）12434x x x ++= 123421x x x x +-+=1234333x x x x --+=-1234234x x x x -++-=实验代码C 语言程序代码如下: #include #include using namespace std;void ColPivot(float *c,int n,float x[]) { int i,j,t,k; float p; for(i=0;i<=n-2;i++){ k=i; for(j=i+1;j<=n-1;j++) if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i)))) k=j; if(k!=j) for(j=i;j<=n;j++){ p=*(c+i*(n+1)+j); *(c+i*(n+1)+j)=*(c+k*(n+1)+j); *(c+k*(n+1)+j)=p; } for(j=i+1;j<=n-1;j++){ p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i)); for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t)); } } for(i=n-1;i>=0;i--){ for(j=n-1;j>=i+1;j--) (*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j)); x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i)); } } int main(){ void ColPivot(float*,int,float[]); int i;float x[4];float c[4][5]={1,1,0,3,4,2,1,-1,1,1,3,-1,-1,3,-3,-1,2,3,-1,4};ColPivot(c[0],4,x);for(i=0;i<=3;i++)printf("[x%d]=%f\n",i,x[i]);return 0;}●运行窗口4、编写用追赶法解三对角线性方程组的程序，并解下列方程组：，其中2）Ax bA10x10=-4 11 -4 11 -4 1. . .. . .1 -4 11 -4b= -27-15…-15●实验代码C语言程序如下：#include#includeusing namespace std;void ColPivot(float *c,int n,float x[]){int i,j,t,k;float p;for(i=0;i<=n-2;i++){k=i;for(j=i+1;j<=n-1;j++)if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i)))) k=j; if(k!=j)for(j=i;j<=n;j++){p=*(c+i*(n+1)+j);*(c+i*(n+1)+j)=*(c+k*(n+1)+j);*(c+k*(n+1)+j)=p;}for(j=i+1;j<=n-1;j++){p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t)); }}for(i=n-1;i>=0;i--){for(j=n-1;j>=i+1;j--)(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));}}int main(){void ColPivot(float*,int,float[]);int i;float x[10];float c[10][11]={-4,1,0,0,0,0,0,0,0,0,-27, 1,-4,1,0,0,0,0,0,0,0,-15,0,1,-4,1,0,0,0,0,0,0,-15,0,0,1,-4,1,0,0,0,0,0,-15,0,0,0,1,-4,1,0,0,0,0,-15,0,0,0,0,1,-4,1,0,0,0,-15,0,0,0,0,0,1,-4,1,0,0,-15,0,0,0,0,0,0,1,-4,1,0,-15,0,0,0,0,0,0,0,1,-4,1,-15,0,0,0,0,0,0,0,0,1,-4,-15};ColPivot(c[0],10,x);for(i=0;i<=9;i++)printf("[x%d]=%f\n",i,x[i]); return 0;}●运行窗口实习题四123●实验代码C语言程序如下：#include#includeusing namespace std;#define N 5void Difference(float x[],float y[],int n){float *f=new float[n+1];int k,i;for(k=1;k<=n;k++){f[0]=y[k];for(i=0;i<k;i++)< p="">f[i+1]=(f[i]-y[i])/(x[k]-x[i]);y[k]=f[k];}delete f;return;}int main(){int i;float a,b,c,varx=0.46,vary=0.55,varz=0.60;float x[N+1]={0.30,0.42,0.50,0.58,0.66,0.72};floaty[N+1]={1.04403,1.08462,1.11803,1.15603,1.19817,1.23223};Difference(x,y,N); a=y[N];b=y[N];c=y[N]; for(i=N-1;i>=0;i--) a=a*(varx-x[i])+y[i]; for(i=N-1;i>=0;i--) b=b*(vary-x[i])+y[i]; for(i=N-1;i>=0;i--) c=c*(varz-x[i])+y[i];printf("Nn(%f)=%f\n",varx,a); printf("Nn(%f)=%f\n",vary,b); printf("Nn(%f)=%f\n",varz,c); return 0;}● 运行窗口实习题六1、用复化梯形公式和复化辛卜生公式计算积分I 1（f ）=?+202x cos 1πdx 。

计算方法上机实验报告上课时间: 2014-2015学年秋学期, 6～14周一. 拉格朗日插值------------------------------------------------------1二. 牛顿插值------------------------------------------------------------3三. 改进欧拉法---------------------------------------------------------5四. 四阶龙格-库塔-----------------------------------------------------7五. 牛顿迭代------------------------------------------------------------9 六．复化Simpson公式------------------------------------------------11七. Romberg算法------------------------------------------------------14八. Seidel迭代法------------------------------------------------------17九. Gauss列主元消去法----------------------------------------------20一.拉格朗日插值1.程序代码#include<iostream.h>void Lagrange(){int i=0;double a[10],b[10],L,L1,L2,L3,L4,x;cout<<"x=";for(i=0;i<4;i++){cin>>a[i];}cout<<"y=";for(i=0;i<4;i++){cin>>b[i];}cout<<"x=";cin>>x;L1=(x-a[1])*(x-a[2])*(x-a[3])*b[0]/(a[0]-a[1])/(a[0]-a[2])/(a[0]-a[3]);L2=(x-a[0])*(x-a[2])*(x-a[3])*b[1]/(a[1]-a[0])/(a[1]-a[2])/(a[1]-a[3]);L3=(x-a[0])*(x-a[1])*(x-a[3])*b[2]/(a[2]-a[0])/(a[2]-a[1])/(a[2]-a[3]);L4=(x-a[0])*(x-a[1])*(x-a[2])*b[3]/(a[3]-a[0])/(a[3]-a[1])/(a[3]-a[2]);L=L1+L2+L3+L4;cout<<"L="<<L;}void main(){Lagrange();cout<<endl;}2.例子3.运行结果二. 牛顿插值1.程序代码#include <iostream.h>#include<conio.h>void main(){int n,i,j;double A[50][50],*x,*y;cout<<"请输入插值节点数: ";cin>>n;x=new double[n];y=new double[n];cout<<"请输入这"<<n<<"个插值节点（xi,yi）: "<<endl;for(i=0;i<=n-1;i++)cin>>x[i]>>y[i];double K=1,xx,N=0,P;for(i=0;i<=n-1;i++){A[i][0]=x[i];A[i][1]=y[i];}for(j=2;j<=n;j++)for(i=1;i<=n-1;i++){A[i][j]=(A[i][j-1]-A[i-1][j-1])/(A[i][0]-A[i-j+1][0]);}for(i=0;i<=n-1;i++)cout<<"输出第"<<i+1<<"阶差商为: "<<A[i][i+1]<<endl;cout<<"请输入预求值x=";cin>>xx;for(i=0;i<n-1;i++){K*=xx-x[i];N+=A[i+1][i+2]*K;P=A[0][1]+N;}cout<<"Newton插值结果为: y="<<P<<endl;getch();}2.例子3.运行结果三. 改进欧拉法1.程序代码#include<iostream.h>#include<conio.h>double fun(double x,double y){return(-0.9*y/(1+2*x));}void main(){double a,b,*y,h,*x,yp,yc;int n,k;cout<<"常微分方程为y'=-0.9*y/(1+2*x)"<<endl;cout<<"其中0<=x<=1"<<endl;cout<<"初值为y(0)=1"<<endl;cout<<"请输入计算区间(a,b): ";cin>>a>>b;cout<<"请输入步长h: ";cin>>h;cout<<"请输入计算次数: ";cin>>n;y=new double[n];x=new double[n];cout<<"请输入初值y(0)=";cin>>y[0];x[0]=a;for(k=0;k<=n;k++){yp=y[k]+h*fun(x[k],y[k]);yc=y[k]+h*fun(x[k]+h,yp);y[k+1]=0.5*(yp+yc);x[k+1]=x[k]+h;}cout<<"迭代结果为: "<<endl;for(k=0;k<=n;k++)cout<<"x"<<k<<"="<<x[k]<<'\t'<<"y"<<k<<"="<<y[k]<<endl;getch();}2.例子3.运行结果四. 四阶龙格-库塔1.程序代码#include<iostream.h>#include<conio.h>double fun(double x,double y){return(x-y);}void main(){double a,b,*y,h,x,k1,k2,k3,k4;int n,k;cout<<"常微分方程为y'=x-y"<<endl;cout<<"其中0<=x<=1"<<endl;cout<<"初值为y(0)=0"<<endl;cout<<"请输入计算区间(a,b): ";cin>>a>>b;cout<<"请输入步长h: ";cin>>h;cout<<"请输入计算次数: ";cin>>n;y=new double[n];cout<<"请输入初值y(0): ";cin>>y[0];x=a;cout<<"Runge-Kutta迭代法结果为: "<<endl;cout<<"x0="<<x<<'\t'<<"y0="<<y[0]<<endl;for(k=0;k<=n-1;k++){k1=fun(x,y[k]);k2=fun(x+h/2,y[k]+k1*h/2);k3=fun(x+h/2,y[k]+k2*h/2);k4=fun(x+h,y[k]+k3*h);y[k+1]=y[k]+(h/6)*(k1+2*(k2+k3)+k4);cout<<"x"<<k+1<<"="<<x+h<<'\t'<<"y"<<k+1<<"="<<y[k+1]<<endl;x=x+h;}getch();}2.例子3.运行结果五. 牛顿迭代法1.程序代码(C++代码)#include<iostream>#include<cmath>using namespace std;double newtondiedai(double a,double b,double c,double d,double x); int main(){double a,b,c,d;double x=1.5;cout<<"请依次输入方程四个系数: ";cin>>a>>b>>c>>d;x=newtondiedai(a,b,c,d,x);cout<<x<<endl;return 0;}double newtondiedai(double a,double b,double c,double d,double x) {while(abs(a*x*x*x+b*x*x+c*x+d)>0.000001){x=x-(a*x*x*x+b*x*x+c*x+d)/(3*a*x*x+2*b*x+c);}return x;}2.例子3.运行结果六. 复化Simpson公式1.程序代码（C++代码）#include<iostream.h>#include<math.h>double function1(double x)//被积函数{double s;s=x/(4+x*x);return s;}double function2(double x)//被积函数{double s;s=sqrt(x);return s;}double ReiterationOfSimpson(double a,double b,double n,double f(double x))//复化Simpson公式{double h,fa,fb,xk,xj;h=(b-a)/n;fa=f(a);fb=f(b);double s1=0.0;double s2=0.0;for(int k=1;k<n;k++){xk=a+k*h;s1=s1+f(xk);}for(int j=0;j<n;j++){xj=a+(j+0.5)*h;s2=s2+f(xj);}double sn;//和sn=h/6*(fa+fb+2*s1+4*s2);//复化Simpson公式return sn;}main(){double a,b,Result,n;cout<<"请输入积分下限:"<<endl;cin>>a;cout<<"请输入积分上限:"<<endl;cin>>b;cout<<"请输入分割区间数n:"<<endl;cin>>n;cout<<"复化Simpson 公式计算结果:";Result=ReiterationOfSimpson(a, b, n,function1); cout<<Result<<endl;}2.例子 （620(n 3)4x I dx x ==+⎰）3.运行结果七. Romberg算法1.程序代码（C++代码）#include<iostream>#include<cmath>using namespace std;#define f(x) (4/(1+x*x))#define epsilon 0.0001#define MAXREPT 10double Romberg(double aa,double bb) { int m,n;double h,x;double s,q;double ep;double *y =new double[MAXREPT]; double p;h=bb-aa;y[0]=h*(f(aa)+f(bb))/2.0;m=1;n=1;ep=epsilon+1.0;while((ep>=epsilon)&&(m<MAXREPT)) { p =0.0;for(int i=0;i<n;i++){ x=aa+(i+0.5)*h;p=p+f(x);}p=(y[0] + h*p)/2.0;s=1.0;for(int k=1;k<=m;k++){ s=4.0*s;q=(s*p-y[k-1])/(s-1.0); y[k-1]=p;p=q;}p=fabs(q-y[m-1]);m=m+1;y[m-1]=q;n=n+n;h=h/2.0;}return (q);}int main(){double a,b;cout<<"Romberg积分,请输入积分范围a,b:"<<endl; cin>>a>>b;cout<<"积分结果:"<<Romberg(a,b)<<endl;system("pause");return 0;}2.例子（1241I dxx=+⎰）3.运行结果八. Seidel迭代法1.程序代码（C++代码）# include <math.h># include <stdio.h># define max 100# define EPS 1e-6float a[3][3]={{10,-1,-2},{-1,10,-2},{-1,-1,5}}; float b[3]={7.2,8.3,4.2};float x[3]={0,0,0};float y[3];float S(int m){int n;float S=0;float y;for(n=0;n<3;n++){if(m==n){}else{S+=a[m][n]*x[n];}}y=(b[m]-S)/a[m][m];return y;}void main(){int i;int F,T=1,k=1;do{F=0;if(T){T=0;}else{k++;for(i=0;i<3;i++){x[i]=y[i];}}for(i=0;i<3;i++){y[i]=S(i);}for(i=0;i<3;i++){if(fabs(x[i]-y[i])>EPS){F=1;}}printf("%d\n",k);}while(((F==1)&&k!=max));printf("迭代次数:%d\n",k);for(i=0;i<3;i++){printf("x[%d]=%1.6f\n",i,y[i]); }}2.例子（1012237.21102238.31253 4.2x x xx x xx x x--=⎧⎪-+-=⎨⎪--+=⎩）3.运行结果九. Gauss列主元消去法1.程序代码（C++代码）#include <stdio.h>#include<conio.h>#include <math.h>#define max_dimension 20int n;static float a[max_dimension][max_dimension]; static float b[max_dimension];static float x[max_dimension];void main(){int i;int j;int d;int row;float temp;float known_items;float l[max_dimension][max_dimension];printf("请输入阶数:");scanf("%d",&n);printf("\n");printf("请输入系数矩阵的值: ");printf("\n");for(i=0; i<n; i++){ printf("输入第%d行的值:",i+1);for (j=0; j<n; j++){scanf("%f",&a[i][j]);}printf("\n");}printf("请输入常数项的值: ");for(i=0; i<n; i++)scanf("%f",&b[i]);printf("\n");for(d=0; d<n-1; d++){row=d;for(i=d+1; i<n; i++){if(fabs(a[i][d])>fabs(a[row][d]))row=i;}if(row!=d){for(j=d; j<n; j++){temp=a[row][j];a[row][j]=a[d][j];a[d][j]=temp;}temp=b[row];b[row]=b[d];b[d]=temp;}for(i=d+1; i<n; i++){l[i][d]=-a[i][d]/a[d][d];for (j=d; j<n; j++){a[i][j]=a[i][j]+a[d][j]*l[i][d];}b[i]=b[i]+b[d]*l[i][d];}}printf("\n");for (i=n-1; i>-1; i--){known_items=0;for(j=1; j<n-i; j++){known_items=known_items+a[i][i+j]*x[i+j]; }x[i]=(b[i]-known_items)/a[i][i];}printf("方程组的根为:\n\n");for(i=0; i<n; i++)printf("%.5f ",x[i]);printf("\n\n");getch();}2.例子3.运行结果。

东南⼤学计算⽅法实验报告计算⽅法与实习实验报告学院：电⽓⼯程学院指导⽼师：李翠平班级：160093姓名：黄芃菲学号：16009330实习题⼀实验1 拉格朗⽇插值法⼀、⽅法原理n次拉格朗⽇插值多项式为：L n(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+y n l n(x)n=1时，称为线性插值，L1(x)=y0(x-x1)/(x0-x1)+ y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0)n=2时，称为⼆次插值或抛物线插值，精度相对⾼些L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1)⼆、主要思路使⽤线性⽅程组求系数构造插值公式相对复杂，可改⽤构造⽅法来插值。

对节点x i(i=0,1,…,n)中任⼀点x k(0<=k<=n)作⼀n 次多项式l k(x k)，使它在该点上取值为1，⽽在其余点x i(i=0,1,…,k-1,k+1,…,n)上为0，则插值多项式为L n(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+y n l n(x) 上式表明：n 个点x i(i=0,1,…,k-1,k+1,…,n)都是l k(x)的零点。

可求得l k三．计算⽅法及过程：1.输⼊节点的个数n2.输⼊各个节点的横纵坐标3.输⼊插值点4.调⽤函数，返回z函数语句与形参说明程序源代码如下：形参与函数类型参数意义int n 节点的个数double x[n]（double *x）存放n个节点的值double y[n]（double *y）存放n个节点相对应的函数值double p 指定插值点的值double fun() 函数返回⼀个双精度实型函数值，即插值点p处的近似函数值#include#includeusing namespace std;#define N 100double fun(double *x,double *y, int n,double p); void main(){int i,n;cout<<"输⼊节点的个数n：";cin>>n;double x[N], y[N],p;cout<<"please input xiangliang x= "<for(i=0;i>x[i];cout<<"please input xiangliang y= "<for(i=0;i>y[i];cout<<"please input LagelangrichazhiJieDian p= "< cin>>p;cout<<"The Answer= "<system("pause") ;}double fun(double x[],double y[], int n,double p) {double z=0,s=1.0;int k=0,i=0;double L[N];while(k{ if(k==0){ for(i=1;iL[0]=s*y[0];k=k+1;}else{s=1.0;for(i=0;i<=k-1;i++)s=s*((p-x[i])/(x[k]-x[i]));for(i=k+1;iL[k]=s*y[k];k++;}}for(i=0;ireturn z;}五．实验分析n=2时，为⼀次插值，即线性插值n=3时，为⼆次插值，即抛物线插值n=1，此时只有⼀个节点，插值点的值就是该节点的函数值n<1时，结果都是返回0的；这⾥做了n=0和n=-7两种情况3常⽤的是线性插值和抛物线插值，显然，抛物线精度相对⾼些n次插值多项式Ln（x）通常是次数为n的多项式，特殊情况可能次数⼩于n.例如：通过三点的⼆次插值多项式L2(x),如果三点共线，则y=L2(x)就是⼀条直线，⽽不是抛物线，这时L2(x)是⼀次式。