计算方法第二章方程求根上机报告

数值计算⽅法⽅程求根数值计算⽅法实验报告实验内容：⽅程求根实验室：专业班级：学号：姓名：2.⽤MATBAB软件，⽤⼆分法求⽅程f（x）=x^3+4*x^2-10=0在区间[1,2]内根的近似值，为使误差不超过10^-5时所需要的⼆分次数。

function bisection_time(tolerance)a=1;b=2;k=0;while(abs(b-a)>tolerance)c=(a+b)/2;fa=a^3+4*a^2-10;fb=b^3+4*b^2-10;fc=c^3+4*c^2-10;if((fa==0)|(fc==0))disp(k);elseif(fa*fc<0)b=c;k=k+1;elseif(fc*fb<0)a=c;k=k+1;elseif(fb==0)disp(k);endendsoluntion=(a+b)/2;disp(soluntion);disp(k);运⾏结果1.36523176.取x0=1.5，⽤⽜顿迭代法求f（x）=x^3+4*x^2-10=0的跟的近似值function new(tolerance)x0=1.5;k=0;a=x0^3+4*x0^2-10;b=3*x0^2+8*x0;x1=x0-a/b;while(abs(x0-x1)>tolerance)x0=x1;k=k+1;a=x0^3+4*x0^2-10;b=3*x0^2+8*x0;x1=x0-a/b;enddisp(x1);disp(k);运⾏结果1.3652338.弦割法求⽅程f（x）=x^3-3*x^2-x+9=0在区间[-2，-1]内的⼀个实根近似值Xk，使|f(x) |<=10^-5. function xuange(k)x0=-2;x1=-1;t=0;a=x1^3-3*x1^2-x1+9;b=x0^3-3*x0^2-x0+9;x2=x1-a*(x1-x0)/(a-b);while(abs(x1-x0)>k)x0=x1;x1=x2;a=x1^3-3*x1^2-x1+9;b=x0^3-3*x0^2-x0+9;x2=x1-a*(x1-x0)/(a-b);t=t+1;enddisp(x1);disp(t)运⾏结果-1.52510269.⽤艾特肯算法求⽅程f （x ）=x^3+4*x^2+10=0在区间[1,2]内的根的近似值（取X0=1.5，g （x ）=410x ，精确到|Xk+1-Xk|<=10^-5，并与第2,3,6题的相应结果进⾏⽐较。

《数值分析》实验报告实验一方程求根一、实验目的:掌握二分法、Newton法、不动点迭代法、弦截法求方程的根的各种计算方法、并实施程序调试和运行，学习应用这些算法于实际问题。

二、实验内容:二分法、Newton法、不动点迭代法、弦截法求方程的根、程序的调试和运行，给出实例的计算结果。

观察初值对收敛性的影响。

三、实验步骤：①、二分法：定义：对于区间[a，b]上连续不断且f（a）·f（b）<0的函数y=f（x），通过不断地把函数f（x）的零点所在的区间一分为二，使区间的两个端点逐步逼近零点，进而得到零点近似值的方法叫二分法。

实现方法：首先我们设一方程400*(x^4)-300*(x^3)+200*(x^2)-10*x-1=0，并求其在区间[0.1,1]上的根，误差限为e=10^-4。

PS：本方法应用的软件为matlab。

disp('二分法')a=0.1;b=1;tol=0.0001;n0=100;fa=400*(a.^4)-300*(a.^3)+200*(a.^2)-10*a-1;for i=1:n0 p=(a+b)/2;fp=400*(p.^4)-300*(p.^3)+200*(p.^2)-10*p-1;if fp==0||(abs((b-a)/2)<tol)disp('用二分法求得方程的根p=')disp(p)disp('二分迭代次数为：')disp(i)break;end;if fa*fp>0 a=p;else b=p;end;end;if i==n0&&~(fp==0||(abs((b-a)/2)<tol)) disp(n0) disp('次二分迭代后没有求出方程的根')end;程序调试：运行结果：用二分法求得方程的根p=0.1108二分迭代次数为：14②Newton法定义：取定初值x0，找到函数对应的点，然后通过该点作函数切线，交x轴，得到新的横坐标值，然后找函数对应的点，做切线，得到新的横坐标值,重复上述步骤，多次迭代，直到收敛到需要的精度。

（1）对原方程变形为x=arctanX+π,令f(x)=x-arctanx-π,f ’(x)=221x x+ （2）选定初值x 0=4.5,构造牛顿迭代公式2211arctan x x x x x x n n +---=+π（3）在VC++6.0中编写程序如下并运行： #include<stdio.h> #include<math.h> #define ESP 1e-6 #define PI 3.1415926 void main() {double x0,x1,f,f1,a,b; printf("INPUT x0:"); scanf("%lf",&x0); a=floor(x0/PI); do {f=x0-atan(x0)-a*PI; f1=x0*x0/(1+x0*x0); x1=x0-f/f1; b=x0; x0=x1; }while(fabs(x0-b)<=ESP); printf("%f\n",x0); }2实验步骤1.1a)程序编译出错：b)找到出错位置：c)修改为“x1=0.2-0.1*exp(x0)”后调试运行为：调试成功，运行程序得出结果。

1.2a)编译运行出错：b)找到出错位置c)修改为“f=x0-atan(x0)-a*PI;”，调试运行。

运行得出结果。

3实验结论（数据及分析结果）（1）对于方程e x+10x-2=0，输入初值为0时结果如下：因为x∈(0.1),φ’(x)=-0.1e x<1,所以该迭代格式收敛。

（2）对于方程x=tgx，输入初值4.5时运行结果如下:当输入x为100时：牛顿迭代是局部收敛的，故迭代在方程的根的附近是收敛的，所以初值的选择对牛顿迭代的收敛性有影响，若初值选在根的附近则迭代收敛，若初值选择离根远则发散。

4实验小结（收获体会）通过这次实验，基本掌握了利用C语言解决数值计算中的方程求根问题，从最初的编写算法到调试再到得出结果，虽然有困难但是通过翻阅资料，查工具书等，最终顺利完成了任务，同时也加深了对于迭代法的认识。

实验一方程求根一、实验目的用各种方法求任意实函数方程f(x)=0在自变量区间[a,b]上，或某一点附近的实根。

并比较方法的优劣。

二、实验方法（1）二分法对方程f(x)=0在[a,b]内求根。

将所给区间等分，在分点x=(b-a)/2判断是否f(x)=0，若是，则有根x=(b-a)/2.否则,继续判断是否f(a)·f(x)<0,若是，则令b=x,否则令a=x。

重复此过程直至求出方程f(x)=0在[a,b]中的近似根为止。

（2）迭代法将方程f(x)=0等价变换为x=h(x)形式，并建立相应的迭代公式Xk+1=h(Xk)。

（3）牛顿法若已知方程f(X)=0的一个近似根X0,则函数f(X)在点X0附近可用一阶泰勒多项式P1= f (X0) + f’ (X0) (X-X0)来近似，因此方程f(X)=0可近似表示为f(X0)+ f’ (X0) (X-X0)=0.设f’ (X0)≠0，则X= X0- f (X0)/ f’ (X0)，取X作为原方程新的近似根X1，然后将X1作为X0带入上式，迭代公示为：X k+1=X k - f (X k)/ f’ (X k)。

三、实验内容在区间[0,1]上用二分法求方程的近似根，要求误差不超过0.5×10^3。

取初值X0=0，用迭代公式X k+1=(2-e^k)/10,(k=0,1,2,…)求方程e^x+10x-2=0的近似根。

要求误差不超过0.5×10^3。

取初值X0=0，用牛顿迭代法求方程e^x+10x-2=0的近似根。

要求误差不超过0.5×10^3。

四、实验程序1.二分法function x=agui_bisect(fname,a,b,e)fa=feval(fname,a);fb=feval(fname,b);if fa*fb>0 error('两端函数值为同号');endk=0x=(a+b)/2while(b-a)>(2*e)fx=feval(fname,x);if fa*fx<0b=x;fb=fx;elsea=x;fa=fx;endk=k+1x=(a+b)/2end2.迭代法function x=agui_iterative(fname,x0,e)N=100;x=x0;x0=x+2*e;k=0;while abs(x0-x)>e & k<Nk=k+1x0=x;x=feval(fname,x0);disp(x)endif k==N warning('已达最大迭代次数');end3.牛顿法function x=agui_newton(fname,dfname,x0,e)N=100;x=x0;x0=x+2*e;k=0;while abs(x0-x)>e&k<Nk=k+1x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);disp(x)endif k==N warning('已达最大迭代次数');end五、实验结果1.二分法2.迭代法3.牛顿法六、结果分析二分法要循环10次，迭代法要迭代4次，牛顿法要迭代3次才能达到精度为0.5×1^-3的要求，由此可知：二分法方法简单，编程容易，且对函数f（x）的性质要求不高，但其收缩速度较慢，计算量大，因此常被用于精度不高的近似根，或为迭代法求初值。

实验四 方程求根实验一. 实验目的（1）深入理解方程求根的迭代法的设计思想，学会利用校正技术和松弛技术解决某些实际的非线性方程问题，比较这些方法解题的不同之处。

（2）熟悉Matlab 编程环境，利用Matlab 解决具体的方程求根问题。

二. 实验要求用Matlab 软件实现根的二分搜索、迭代法、Newton 法、快速弦截法和弦截法，并用实例在计算机上计算。

三. 实验内容1. 实验题目（1）早在1225年，古代人曾求解方程020102)(23=-++=x x x x f 并给出了高精度的实根368808107.1*=x ，试用Newton 法和弦截法进行验证，要求精度610-=ε，并绘制方程的图形。

答：A.Newton 法：a ．编写文件Newton.m 、func4.m 内容如下所示：b.运行，如下所示A为矩阵，由上面可知，对于初值为5，运行7次即可得到所需的精度，验证结果为古人给出的解释正确的；c.作图，编写下面的文件photo1.m.然后运行即可：注意下面中的x矩阵即为刚才计算出来的x系列，k为迭代的次数：a．编写文件Chord.m内容如下所示：b．运行结果如下所示：由上表可知，在精度为10^-6时有7位有效数字，古人的结果还是正确的c．作图，在上面运行后，即运行newton法时写的photo1.m文件即可出现图像：可以看到图中两条曲线基本重合； （2）取5.00=x ，用迭代法求方程x e x -=的根，然后用Aitken 方法加速，要求精度为结果有4为有效数字。

答：a. 编写文件func7.m 和Aiken.m ，内容如下所示：b ．运行：具有四位有效数字 （3）用快速弦截法求解方程01)(=-=x xe x f ，要求精度为610-=ε，取6.05.010==x x ，作为开始值，并绘制1)(-=x xe x f 的图形。

答：对照可知，书本后面的程序已经正确，运行即可：下面为快速弦截法的主程序文件：函数文件如下：运行如下：作图，编写下面的文件：运行该文件就可以y=x*exp(x)-1函数和插值函数的图：可以看到两条直线基本重合在一起了，扩大图片可以看到两条直线是不重合的：2. 设计思想要求针对上述题目，详细分析每种算法的设计思想。

方程求根§2.0 引言§2.1 二分法§2.2 简单迭代法§2.3 牛顿（Newton）法§2.4 其它求根方法（迭代过程的加速方法）§2.5 作业讲评2.0 引 言非线性科学是当今科学发展的一个重要研究方向，非线性方程的求根也成为其中一个重要内容。

一般而言，非线性方程的求根非常复杂。

在实际应用中有许多非线性方程的例子，例如（1）在光的衍射理论（the theory of diffraction of light)中,需要求x-tanx=0的根（2）在行星轨道（ planetary orbits ）的计算中,对任意的a 和b ，需要求x-asinx=b 的根（3）在数学中，需要求n 次多项式-1-110 ... 0n n n n a x a x a x a ++++=的根。

非线性方程的一般形式 ()0f x = 这里()f x 是单变量x 的函数，它可以是代数多项式-1-110() ... nn n n f x a x a x a x a =++++ (0n a ≠)也可以是超越函数，即不能表示为上述形式的函数。

满足方程 ()0f x = 的x 值通常叫做方程的根或解，也叫函数()0f x =的零点。

2.1 二分法(Bisection Method)1 概念:二分法也称对分区间法、对分法等，是最简单的求根方法，属于区间法求根类型。

在用近似方法时，需要知道方程的根所在区间。

若区间[a,b]含有方程f(x)=0的根，则称[a,b]为f(x)=0的有根区间；若区间[a,b]仅含方程f(x)= 0的一个根，则称[a,b]为f(x)= 0的一个单根区间。

2．基本思想根的存在定理(零点定理)：f(x)为[a,b]上的连续函数，若f(a)·f(b)<0，则[a,b]中至少有一个实根。

如果f(x)在[a,b]上还是单调递增或递减的，则f(x)=0仅有一个实根。

实验报告名称方程求根班级：学号：姓名：成绩：1实验目的（1）通过对二分法与牛顿迭代法做编程练习和上机运算，进一步体会二分法和牛顿法的不同。

（2）编写割线迭代法的程序，求非线性方程的解，并与牛顿迭代法作比较。

2 实验内容求方程f(x)=x^3+x^2-3*x-3在1.5附近的根。

3实验步骤二分法算法给定区间（a,b），并设f(a)与f(b)符号相反，取ε为根的容许误差，δ为|f(x）|的容许误差。

①令c=（a+b）/2.②如果（c-a）<ε或|f（c）|<δ,则输出c，结束；否则执行③。

③如果f（a）*f(c)>0,则a=c；否则令b=c，重复①②③。

牛顿迭代法算法给定初始值x0，ε为根的容许误差，η为|f(x)|的容许误差，N为迭代次数的容许值。

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

②计算x1=x0-f(x0)/f(x0)’.③令x0=x1，转向①。

4 程序设计二分法c语言程序设计：#include <stdio.h>#include <math.h>#define eps 5e-6#define delta 1e-6float Bisection(float a,float b,float(*f)(float)){float c, fc,fa=(*f)(a),fb=(*f)(b);int n=1;printf("二分次数\t\tc\t\tf(c)\n");while(1){if(fa*fb>0){printf("不能用二分法求解");break;}c=(a+b)/2,fc=(*f)(c);if(fabs(fc)<delta)break;else if(fa*fc<0){ b=c;fb=fc;}else{a=c;fa=fc;}if(b-a<eps)break;printf("%d\t\t%f\t\t%f\n",n++,c,fc);}return c;}float f(float x){return x*x*x+x*x-3*x-3;}void main(){float a=1,b=2；float x;x=Bisection(a,b,f);printf("\n方程的根是：%f",x);}牛顿法c语言程序设计：#include<stdio.h>#include<math.h>#define N 100#define eps 1e-6#define eta 1e-8float Newton(float(*f)(float),float(*f1)(float),float x0) {float x1,d;int k=0;do{x1=x0-(*f)(x0)/(*f1)(x0);if(k++>N||fabs((*f1)(x1))<eps){printf("\n Newton迭代散发");break;}d=fabs(x1)<1?x1-x0:(x1-x0)/x1;x0=x1;printf("x(%d)=%f\t",k,x0);}while(fabs(d)>eps&&fabs((*f)(x1))>eta);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 x0,y0;printf("请输入迭代值x0\n");scanf("%f",&x0);printf("x(0)=%f\n",x0);y0=Newton(f,f1,x0);printf("方程的根为:%f\n",y0);}5实验结果及分析二分法的输出结果：牛顿法的输出结果：实验分析：上面程序取三个不同初值，得到同样的结果，但迭代次数不同，初值越接近所求的根，迭代次数越少。

实验报告名称班级：学号：姓名：成绩：1实验目的1）通过对二分法与牛顿迭代法作编程练习与上级运算，进一步体会二分法与牛顿迭代法的不同特点。

2）编写割线迭代法的程序，求非线性迭代法的解，并与牛顿迭代法。

2 实验内容用牛顿法和割线法求下列方程的根x^2-e^x=0;x*e^x-1=0;lgx+x-2=0;3实验步骤1)根据二分法和牛顿迭代法，割线法的算法编写相应的求根函数；2)将题中所给参数带入二分法函数，确定大致区间；3)用牛顿迭代法和割线法分别对方程进行求解；3 程序设计牛顿迭代法x0=1.0;N=100;k=0;eps=5e-6;delta=1e-6;while(1)x1=x0-fc1(x0)/fc2(x0);k=k+1;if k>Ndisp('Newmethod failed')breakendif(abs(x1-x0)<delta || abs(fc1(x1))<delta)break;endx0=x1; %²»ÄÜ·ÅÔÚabs(x1-x0)Ç°endfprintf('%f',x0)fprintf(' %f ', abs(fc1(x1)) )割线法function cutline(x0,x1)N=100;k=0;delta=5e-8;while(1)while(abs(x1-x0)>=delta)c=x1;x1=cutnext(x0,x1);x0=c; %x0 x1µÝÍƵõ½x1 x2 ÈÔÈ»±£´æÔÚx0 x1 endk=k+1;if k>Ndisp('Cutline method failed')break;endif(abs(x1-x0)<delta || abs(fc1(x1))<delta)break;endendfprintf('%.10f\n',x1);function y=cutnext(a,b)y=b-fc(b)/(fc(b)-fc(a))*(b-a);1）原函数function fc1=fc1(x)fc1=x^2-exp(x);end导函数function fc2=fc2(x)fc2=2*x-exp(x);end2）原函数function fc1=fc1(x)fc1=x*exp(x)-1;end导数function fc2=fc2(x)fc2=(x+1)*exp(x);end3）原函数function fc1=fc1(x)fc1=log10(x)+x-2;end导函数function fc2=fc2(x)fc2=1/x/log(10)+1;end4实验结果及分析1）牛顿法结果-0.7034722378割线法结果-0.70346742252）牛顿法结果0.5671435302割线法结果0.56714329043）牛顿法结果1.7553985566割线法结果1.7555794993牛顿迭代法由于设置delta=1e-6，所以算出的误差e<1.0*10^-5; 割线法由于设置delta=5e-8，所以误差e<1.0*10^-7;5总结编程时由于将迭代的代码x0=x1放在if(abs(x1-x0)<delta || abs(fc1(x1))<delta)break;end之前导致程序没有执行就跳出，通过Debug发现了问题，将x0=x1；放到了循环体内部的最后一行，程序得以成功的运行。

简单迭代法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(12*x+sin(x)-1,1.0/3)void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果牛顿迭代法一#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)-sin(x)-12*x+1)/(3*pow(x,2)-cos(x)-12) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果埃特金加速法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) (pow(x,3)-sin(x)+1)/12void main(){int i;double x1=x0,x2=x0;double z,y;printf("k\tx1\tx2\txk\n");for(i=0;i<MAXREPT;i++){if(i==0)printf("%d\t\t\t%g\n",i,x2);elseprintf("%d\t%g\t%g\t%g\n",i,y,z,x2);y=G(x1);z=G(y);x2=z-((z-y)*(z-y))/(z-2*y+x1);if (fabs(x2-x1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x2,i);return;}x1=x2;}printf("AFTER %d repeate,no solved.\n",MAXREPT);} 结果牛顿迭代法二#include<stdio.h>#include<math.h>#define x0 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)+pow(x,2)-3*x-3)/(3*pow(x,2)+2*x-3) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果弦截法#include<stdio.h>#include<math.h>#define x0 0#define x1 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(x,3)+pow(x,2)-3*x-3void main(){int i;double x_k=x0,x_k1=x1,x_k2=0;double y,z;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k2);y=G(x_k);z=G(x_k1);x_k2=x_k1-(z*(x_k1-x_k))/(z-y);if (fabs(x_k2-x_k1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k2,i);return;}x_k=x_k1;x_k1=x_k2;}printf("AFTER %d repeate,no solved.\n",MAXREPT); } 结果高斯顺序消元法#include<stdio.h>#include<math.h>#define N 4static double aa[N][N+1]={{2,4,0,1,1},{3,8,2,2,3},{1,3,3,0,6},{2,5,2,2,3}}; int gauss(double a[][N+2],double x[]);void putout(double a[][N+2]);void main(){int i,j,det;double a[N+1][N+2],x[N+1];for(i=1;i<=N;i++)for(j=1;j<=N+1;j++)a[i][j]=aa[i-1][j-1];det=gauss(a,x);if(det!=0)for(i=1;i<=N;i++)printf(" x[%d]=%g",i,x[i]);printf("\n");}int gauss(double a[][N+2],double x[]){int i,j,k;double c;putout(a);for(k=1;k<=N-1;k++){ if(fabs(a[k][k])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(i=k+1;i<=N;i++){c=a[i][k]/a[k][k];for(j=k;j<=N+1;j++){a[i][j]=a[i][j]-c*a[k][j];}}putout(a);}if(fabs(a[N][N])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(k=N;k>=1;k--){x[k]=a[k][N+1];for(j=k+1;j<=N;j++){x[k]=x[k]-a[k][j]*x[j];}x[k]=x[k]/a[k][k];}return(1);}void putout(double a[][N+2]){for(int i=1;i<=N;i++){for(int j=1;j<=N+1;j++)printf("%-15g",a[i][j]);printf("\n");}printf("\n");}结果雅克比迭代法#include<stdio.h>#include<math.h>#define N 5#define EPS 0.5e-4static double aa[N][N]={{4,-1,0,-1,0},{-1,4,-1,0,-1},{0,-1,4,-1,0},{-1,0,-1,4,-1},{0,-1,0,-1,4}}; static double bb[N]={2,1,2,1,2};void main(){int i,j,k,NO;double a[N+1][N+1],b[N+1],x[N+1],y[N+1];double d,sum,max;for(i=1;i<=N;i++){for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];}printf("\n 请输入最大迭代次数（尽量取大值！）NO:");scanf("%d",&NO);printf("\n");for(i=1;i<=N;i++)x[i]=0;k=0;printf(" k",' ');for(i=1;i<=N;i++)printf("%8cx[%d]",' ',i);printf("\n 0");for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");do{for(i=1;i<=N;i++){sum=0.0;for(j=1;j<=N;j++)if(j!=i) sum=sum+a[i][j]*x[j];y[i]=(-sum+b[i])/a[i][i];}max=0.0;for(i=0;i<=N;i++){d=fabs(y[i]-x[i]);if(max<d) max=d;x[i]=y[i];}printf("%6d",k+1);for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");k++;}while((max>=EPS)&&(k<NO));printf("\nk=%d\n",k);if(k>=NO) printf("\nfail!\n");elsefor(i=1;i<=N;i++)printf("x[%d]=%g\t",i,x[i]);}结果拉格朗日插值多项式#include<stdio.h>#include<math.h>#define N 4doublex[N]={0.56160,0.56280,0.56401,0.56521},y[N]={0.82741,0.82659,0.82577,0.82495}; void main(){double x=0.5635;double L(double xx);double lagBasis(int k,double xx);void output();output();printf("\n近似值L(%g)=%g\n",x,L(x));}double lagBasis(int k,double xx){double lb=1;int i;for(i=0;i<N;i++)if(i!=k) lb*=(xx-x[i])/(x[k]-x[i]);return lb;}double L(double xx){double s=0;int i;for(i=0;i<=N;i++)s+=lagBasis(i,xx)*y[i];return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<N;i++)printf("\t%g",x[i]);printf("\nyi:");for(i=0;i<N;i++)printf("\t%g",y[i]);}结果牛顿插值多项式#include <math.h>#include <stdio.h>int a;#define M 4double x[M+1]={0.4,0.55,0.65,0.8,0.9},y[M+1]={0.41075,0.57815,0.69675,0.88811,1.02652}; void main(){double x;printf("输入x=");scanf("%lf",&x);printf("次数：");scanf("%d",&a);double N(double xx,int a);void output();output();printf("\n%d次牛顿插值多项式N(%g)=%g\n",a,x,N(x,a));}double N(double xx,int a){double s=y[0],d=1;int i,j;double df[M+1][M+1];for(i=0;i<=M;i++)df[i][0]=y[i];for(j=1;j<=a;j++)for(i=j;i<=a;i++)df[i][j]=(df[i][j-1]-df[i-1][j-1])/(x[i]-x[i-j]);printf("\nx\tf(x)\t");for(j=1;j<=a;j++) printf("%5d阶",j);for(i=0;i<=a;i++){printf("\n%g\t%g",x[i],y[i]);for(j=1;j<=i;j++)printf("\t%7.5g",df[i][j]);}for(i=1;i<=a;i++){d*=(xx-x[i-1]);s+=df[i][i]*d;}return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<=M;i++)printf("\t%7g",x[i]);printf("\nyi:");for(i=0;i<=M;i++)printf("\t%7g",y[i]);}结果复合梯形公式#include<stdio.h>#include<math.h>#define f(x) 1/(x*x+1)#define Pi 3.1415926void main(){double a=0,b=1;double T,h,x;int n,i;printf("please input n:");scanf("%d",&n);h=(b-a)/n;x=a;T=0;for(i=1;i<n;i++){x+=h;T+=f(x);}T=(f(a)+2*T+f(b))*h/2;printf("T(%d)=%g\n",n,T);printf("The exact value is %g\n",Pi/4);}复合辛普森公式#include<stdio.h>#include<math.h>#define f(x) 1/(1+x*x)#define Pi 3.1415926void main(){double a=0,b=1;double S,h,x;int n,i;printf("please input Even n:");scanf("%d",&n);h=(b-a)/n;x=a; S=0;for(i=1;i<n;i++){x+=h;if(i%2==0) S+=2*f(x);else S+=4*f(x);}S=(f(a)+S+f(b))*h/3;printf("S(%d)=%g\n",n,S);printf("The exact value is %g\n",Pi/4);}龙贝格公式加速#include<stdio.h>#include<math.h>#define f(x) sin(x)/(1+x)#define M 3void main(){double a=0,b=1;double Long(double a,double b);printf("近似值I=%g\n",Long(a,b));}double Long(double a,double b){int n=1,i=1,j=1;double T[M+1][M+1],h,x,sum;h=b-a;T[0][0]=(f(a)+f(b))/2;for(j=1;j<=3;j++){x=a;h/=2;n*=2;sum=0;for(i=1;i<=n;i+=2){x=a+i*h;sum+=f(x);}T[j][0]=T[j-1][0]/2+h*sum;}for(i=1;i<=M;i++)for(j=1;j<=i;j++){T[i][j]=(pow(4,j)*T[i][j-1]-T[i-1][j-1])/(pow(4,j)-1);}printf("k\tT0\tT1\tT2\tT3\n");for(i=0;i<=M;i++){printf("%d",i);for(j=0;j<=i;j++)printf(" %g",T[i][j]);printf("\n");}return T[M][M];}。

第二章方程求根2.1 方程求根与二分法2.1.1 引言单个变量的方程(2.1.1)求根是数值计算经常遇到的问题.当f(x)为一般连续函数时，称式(2.1.1)为超越方程，如果f为多项式若,为n次多项式，此时方程(2.1.1)称为代数(或多项式)方程.如果x*(实数或复数)使，则称x*为方程(2.1.1)的根，若，m≥1，且，当m＞1时，称x*为方程(2.1.1)的m重根或称x*是f的m重零点.若x*是f的m重零点，且g充分光滑，则。

当f为式(2.1.2)表示的代数多项式时，根据代数基本定理可知方程(2.1.1)有n个根(含复根，m重根为m个根)，对n=2的代数方程它的根可由公式表示为而当n=3,4时方程的根虽可用公式表示，但表达式太复杂，一般不用，当n≥5已没有直接用公式表达的求根算法.因此对n≥3的代数方程求根方法与一般超越方程(2.1.1)一样都采用迭代方法求根，设(表示f在区间上连续)，若有f()f(b)＜0,则f(x)=0在区间上至少有一个实根，称为有根区间，通常可用逐次搜索法求得方程(2.1.1)的有根区间.例2.1 求方程的有根区间.解根据有根区间定义，对方程的根进行搜索计算，结果如下表：方程的三个有根区间为［1,2］,［3,4］,［5,6］.讲解：非线性方程f(x)=0求根，包括求超越方程和代数方程的根x*,方程的根也是f(x)的零点，即f(x*)=0,x*可以是实根也可以是复根，本章以求实根为主。

求实根首先要确定根x*所在区间，称为有根区间。

根据连续函数性质，若f(x)在上连续，当f()f(b)<0时，为有根区间,为找到方程f(x)=0的有根区间,可用逐次搜索法,也就是在x的不同点上计算f(x),观察f(x)的符号,如例2.1表中所示,只要在相邻两点f反号,则得到有根区间,本例得到三个有根区间,分别为[1,2][3,4][5,6].2.1.2 二分法设,且为有根区间，取中点，将它分为两半，检查与是否同号，若是，说明根x*仍在右侧，取,否则取,得到新的有根区间长度仅为的一半(见图2-1).重复以上过程，即取，将再分半，确定根在的哪一侧，得到新区间，其长度为的一半，从而可得一系列有根区间其中每一个区间长度都是前一个区间长度的一半，因此，的长度为且,即为方程(2.1.1)的根x*的一个足够精确的近似根，且误差为(2.1.3)以上过程称为解方程的二分法.它计算简单且收敛性有保证，但收敛较慢，通常可用于求迭代法的一个足够好的初始近似.例2.2求方程在区间［1.0,1.5］内的一个实根，要求准确到小数点后第二位.解这里=1.0,b=1.5，而f()＜0,f(b)＞0，故在［1.0,1.5］中有根，由式(2.1.3)知即,当n=6时得,已达到,精度要求，各次计算结果见表2-1.故为方程的近似根，误差不超过0.005.讲解：从有根区间出发用二分法将它逐次分半，得到的近似根,其误差为(2.1.3)式．它是求实根近似的有效方法，但由于收敛慢，只用作迭代法的初始近似．2.2 迭代法及其收敛性2.2.1 不动点迭代法求方程(2.1.1)的根时将方程改写为等价形式(2.2.1)若x*满足，则称x*为φ的一个不动点，x*也是方程(2.1.1)的一个根，如果φ连续，可构造迭代法(2.2.2)称为不动点迭代法，φ称为迭代函数.若给定初始近似，由式(2.2.2)逐次迭代得到序列，如果，则由式(2.2.2)两端取极限，得即x*为φ的不动点.从几何图象(参见图2-2)可知，（参见图2-2）可知，φ的不动点就是直线与曲线的交点，横坐标x*即为不动点，从它的某个初始近似出发，在曲线确定一点，引平行x轴直线，与直线交于点，其横坐标即为,由式(2.2.2)逐次求得即为如图2-2所示点的横坐标.例2.3求方程在附近的根.解若将方程改写为,构造迭代法(2.2.3) 由可知，显然不收敛.若将方程改为，构造迭代法(2.2.4)计算结果见表2-2.表 2－2从结果看它是收敛的，且在6位有效数字时即为根x*的近似值.例题表明构造迭代法(2.2.2)，必须使迭代序列｛｝收敛，才能求得方程(2.2.1)的解x*.下面给出解的存在唯一性和迭代收敛性定理.定理2.1设，如果对有≤φ(x)≤b,且存在常数L∈(0,1)使(2.2.5)则φ在区间上存在唯一不动点x*，且由式(2.2.2)生成的迭代序列｛｝对任何收敛于x*，并有误差估计(2.2.6)证明先证明不动点存在性，记,由定理条件有及,若有一等号成立，则或,即φ有不动点，否则必有，因故必有,使,x*即为φ的不动点.再证明唯一性，设都是φ的不动点，且则由式(2.2.5)得与假设矛盾，这表明，即不动点是唯一的.下面证明由式(2.2.2)生成的迭代序列｛｝收敛于φ的唯一不动点x*，由于φ(x)∈，故,再由式(2.2.5)有因0＜L＜1,故,即再利用式(2.2.5)考察上式中令则得式(2.2.6).定理证毕.推论若(表示φ在上一阶导数连续)，则定理2.1中的条件式(2.2.5)可改为(2.2.7)则定理2.1中结论全部成立.实际上，由微分中值定理可得，对有故式(2.2.5)成立.以后使用时，如果连续都可用式(2.2.7)代替式(2.2.5).根据定理2.1可验证例2.3中迭代(2.2.4)收敛，因迭代函数，在［1,2］中,故迭代序列｛｝收敛.且,,而对迭代法(2.2.3)，因在区间［1,2］中，故迭代(2.2.3)是发散的.讲解：求方程f(x)=0的不动点迭代法首先要构造迭代序列,就是将方程改写为便于迭代的形式(2.2.1)，同一方程可以写成多种形式，但必须采用使迭代（2.2.2）收敛的迭代函数,如例2.3中f(x)可改写成,也可以写成,(这里及,前者不收敛，后者收敛，如图2.2所示，就是收敛的情形，怎样判断迭代法（2.2.2）是否收敛就要根据定理2.1检验在区间上是否满足以下两个条件：（1）对有（2）在上满足条件（2.2.7）。