非线性方程的数值计算方法实验

1、数值计算方法非线性方程的数值计算方法实验一、 实验描述：在科学研究和工程实践中，经常需要求解大量的非线性方程。本实验正是通过计算机的程序设计，使用迭代法、波尔查诺二分法、试值法、牛顿-拉夫森法和割线法，来实现非线性方程的求解。本实验中通过对各种方法的实践运用，可以比较出各种方法的优缺点。并且，通过完成实验，可加深对各种方法的原理的理解，熟悉掌握C语言在这些方法中的运用。二、实验内容：1、 求函数的不动点（尽可能多）近似值，答案精确到小数点后12位；2、 如果在240个月内每月付款300美元，求解满足全部年金A为500000美元的利率I，的近似值（精确到小数点后10位）。3、 利用加速牛顿-拉夫

2、森算法，用其求下列函数M阶根p的近似值。(a)、f(x)=(x-2)5，M=5,p=2，初始值p0=1。(b)、f(x)=sin(x3)，M=3,p=0，初始值p0=1。(c)、f(x)=(x-1)ln(x),M=2,p=1，初始值p0=2。 4、 设投射体的运动方程为： y=f(t)=9600(1-e-t/15)-480t x=r(t)=2400(1-e-t/15)(a)求当撞击地面时经过的时间，精确到小数点后10位。 (b)求水平飞行行程，精确到小数点后10位。三、实验原理： (1)、不动点迭代法：它是一种逐次逼近的方法,即用某个固定公式反复校正根的近似值,使之逐步精确化,最后得到满足精度

3、要求的结果。它利用计算机运算速度快，适合做重复性操作的特点，让计算机对一个函数进行重复执行，在每次执行这个函数时，都从变量的原值推出它的一个新值，直至推出最终答案为止。 迭代法一般可用于寻找不动点，即：存在一个实数P，满足P=g(P)，则称P为函数g(x)的一个不动点。且有定理：若g(x)是一个连续函数，且pnn=0是由不动点迭代生成的序列。如果limnpn=P,则P是g(x)的不动点。所以，不动点的寻找多用迭代法。(2)、波尔查诺二分法：起始区间a,b必须满足f(a)与f(b)的符号相反的条件。由于连续函数y=f(x)的图形无间断，所以它会在零点x=r处跨过x轴，且r在区间内。通过二分法可将

4、区间内的端点逐步逼近零点，直到得到一个任意小的包含零点的间隔。二分法定理：设fC(a,b),且存在数ra,b满足f(r)=0。如果f(a)和f(b)的符号相反，且cnn=0为二分法生成的中点序列，则：r-cnb-a2n+1 其中n=0,1, (1)这样，序列cnn=0收敛到零点x=r即可表示为：limncn=r (2)(3)、试值法：假设一个函数中，有f(a)和f(b)符号相反。二分法使用区间a,b的中点进行下一次迭代。如果找到经过点(a,f(a)和(b,f(b)的割线L与x轴的交点(c,0),则可得到一个更好的近似值。为了寻找值c，定义了线L的斜率m的两种表示方法，一种表示方法为：m=fb-

5、f(a)b-a (3)这里使用了点(a,f(a)和(b,f(b)。另一种表示方法为：m=0-f(b)c-b (4)这里使用了点(c,0)和(b,f(b)。 使式(3)和式(4)的斜率相等，则有：fb-f(a)b-a=0-f(b)c-b (5)为了更容易求解c，可进一步表示为：c=b-fb(b-a)fb-f(a) (6)这样会出现3种可能性： 如果f(a)和f(c)的符号相反，则在a,c内有一个零点。 如果f(c)和f(b)的符号相反，则在c,b内有一个零点。如果f(c)=0，则c是零点。然后，可按二分法的方法进行下一步运算。(4)、牛顿-拉夫森法： 此法根据，牛顿-拉夫森定理：设fC2a,b,

6、且存在数pa,b,满足f(p)=0。如果f(p)0,则存在一个数>0，对任意初始近似值p0p-,p+，使得由如下迭代定义的序列pkk=0收敛到p: pk=g(pk-1)= pk-1f(pk-1)f(pk-1) 其中k=1,2, (7)其中，函数g(x)由如下定义：g(x)=x-f(x)f(x) (8)且被称为牛顿-拉夫森迭代函数。由于f(p)=0,显然g(p)=p。这样，通过寻找函数的不动点，可以实现寻找方程f(x)=0的根的牛顿-拉夫森迭代。附：定理(牛顿-拉夫森迭代的加速收敛)： 设牛顿-拉夫森算法产生的序列线性收敛到M阶根x=p,其中M>1,则牛顿-拉夫森迭代公式： pk=p

7、k-1-M f(pk-1)f(pk-1) (9)(5)、割线法： 割线法包含的公式与试值法的公式一样，只是在关于如何定义每个后续项的逻辑判定上不一样。需要两个靠近点(p,0)的初始点(p0,f(p0)和(p1,f(p1)。 可根据两点迭代法公式，得到一般项: pk+1=g(pk, pk-1)= pk-fpk(pk-pk-1)fpk-f(pk-1) (10)四、结果计算及分析：1、函数g(x)=xx-cos(x)的不动点的迭代(迭代法)： (a)计算结果：(b)结果分析:此题经过matlab图形仿真，可以看出其解的大致范围在0.8,1.2之间。 此结果是在取p0=0.879的情况下得出的，其误差

8、精度取为0.000000000001。从迭代过程的误差收敛速度可以看出不动点迭代的误差收敛较慢，所需迭代次数较多。2、满足条件的利率I的计算(试值法)：(a)、计算结果：(b)、结果分析：此题经过一定的计算分析，将书中所给公式结合题中的条件可得出最终的计算式：3600I1+I12240-1-500000=0,求解其中的I即为利率。又经过一定的计算分析，取其初始计算区间为0.1,0.2，误差精度取为0.000000000001。从其误差收敛的速度可以看出，试值法的误差收敛速度快，所需的迭代次数少。3、利用加速牛顿-拉夫森算法，求函数f(x)=(x-2)5、f(x)=sin(x3)和f(x)=(x

9、-1)ln(x)的M阶根p的近似值。其中每个函数的M、p0和最终结果p已给出。(a)、计算结果：(b)、结果分析： 由于C语言函数库中没有求导功能函数，而此题所要求用的加速牛顿-拉夫森算法的计算公式 pk=pk-1-M f(pk-1)f(pk-1) 中要用到题中所给函数的导函数，所以求出它们的导函数分别为f(x)=5(x-2)2, f(x)=3x2cos(x3)和f(x)=In(x)-1- 1x 。此题的误差精度取为0.00000000001。从输出结果中可以大致判断出，加速牛顿-拉夫森法的收敛速度非常快。4、给出投射物体的运动方程y=f(t)=9600(1-e-t/15)-480t和x=r(

10、t)=2400(1-e-t/15)求物体落地时经过的时间和水平飞行程。(割线法) (a)计算结果：(b)、结果分析：此题用了割线法的思路进行程序设计，由于割线法的公式要求，需要给出t0,t1 两个初始值。经过一定的分析，给出t0=5.0，t1=15.0，误差精度设为0.0000000001。从输出结果中可以看出，割线法的误差收敛速度比较快，所需的迭代次数比较少。五、实验结果分析： 经过以上四个实验的比较分析，结合书上的实例与原理分析，可以得出。不动点迭代法比起其他的迭代法，迭代速度慢、次数多。试值法的迭代次相对少一些，但若对其初始区间估算得太大的话，会导致迭代步数增多，甚至会导致发散。割线法比

11、试值法要快得多，可以达到相对较高的精度，且在其迭代过程中每步只需一次新的函数赋值，其收敛阶一般能达到1.1618。加速牛顿-拉夫森法再求多重根时速度很快，收敛阶能达到2。附件：第一题：#include<stdio.h>#include<math.h>#include<float.h>/为了调用FLT_EPSILON，防止出现分母为零的情况。/#define N 1000/保证有足够多的运算次数能达到比较精确的结果，且便于修改。/#define pre 0.0000000000001/使循环能够在达到某一误差精度后停下来。/double main()doubl

12、e pN,m,errN,relerrN;int k,n;p0=0.879;/给定一个合理的初始值以此来进行迭代。/for(k=0;k<N;k+)m=pk-cos(pk);pk+1=pow(pk,m);/进行迭代的方程，由书上给出。/errk=fabs(pk+1-pk);/对绝对误差的运算。/relerrk=errk/(fabs(pk+1)+FLT_EPSILON);/对相对误差的运算。/n=k;/把k赋给n,使下面输出时能够在合适的时候停下来。/if(errk<pre|relerrk<pre) break;/判断循环在何处达到规定精度，并跳出循环。/for(k=0;k<

13、n;k+)printf("nX%d=%14.12f err%d=%14.12f relerr%d=%14.12f",k,pk,k,errk,k,relerrk);/对迭代点及误差的输出。/getch();return 0;第二题：#include<stdio.h>#include<math.h>#include<float.h>#define N 10000#define pre 0.000000000001/使循环能够在达到某一误差精度后停下来。/double main()double f(double x);double aN,bN,

14、IN,errN,relerrN;int i,n;a0=0.10,b0=0.20;/预估它们的初始区间/for(i=0;i<N;i+)Ii=bi-(f(bi)*(bi-ai)/(f(bi)-f(ai);/试值法公式，有书上给出。/ if(f(ai)*f(Ii)>0)ai+1=Ii,bi+1=bi;/根据试值法，判断区间。/else ai+1=ai,bi+1=Ii;/根据试值法，判断区间。/ erri=fabs(Ii-Ii-1);/对绝对误差的运算。/ relerri=erri/(fabs(Ii)+FLT_EPSILON);/对相对误差的运算。/ n=i;/将i值赋给n,使输出能够在合

15、适的时候停下来。/ if(erri<pre|relerri<pre) break;/判断循环在何处达到规定精度，并跳出循环。/ err0=I0-0;/定义初始误差，便于下面输出。/ relerr0=err0/(fabs(I0)+FLT_EPSILON);/定义初始误差，便于下面输出。/ for(i=0;i<n;i+) printf("n I%d=%12.10f err%d=%12.10f relerr%d=%12.10f",i,Ii,i,erri,i,relerri);/输出利率的计算过程。/ getch(); return 0;double f(doub

16、le x) return(3600/x)*(pow(1+x/12,240)-1)-500000);/此题的计算公式，由书上给出。/第三题：#include<stdio.h>#include<math.h>#include<float.h>/为了调用FLT_EPSILON，防止出现分母为零的情况。/#define pre 0.00000000001/使循环能够在达到某一误差精度后停下来。/#define N 1000/迭代次数，以保证能得到最精确的值。/double main() double fa(double x);/声明fa函数。/ double fb(

17、double x);/声明fb函数。/ double fc(double x);/声明fc函数。/ double aN,bN,cN,errN;/定义结果数组，与误差数组。/ int k,i,j,m; a0=1.0,b0=1.0,c0=2.0;/赋初始值，由书上给出。/ for(k=0;k<N;k+)ak+1=ak-5.0*fa(ak);/加速牛顿-拉夫森公式，见书上P64。下同。/errk=fabs(2.0-ak+2);/对绝对误差的运算。k加2，以保证下面输出时其有足够多的值输出。/i=k;/将循环停止时的k值赋给i，便于在下面输出时有恰当的值使输出停止。下同。/if(errk<

18、pre)break;for(k=0;k<N;k+)bk+1=bk-3.0*fb(bk);errk=fabs(0.0-bk);/对绝对误差的运算。/j=k;if(errk<pre)break;for(k=0;k<N;k+)ck+1=ck-2.0*fc(ck);errk=fabs(1.0-ck);/对绝对误差的运算。/m=k;if(errk<pre)break; if(i>j) j=i;if(m>j) m=j;/对i,m,j进行比较，并求出最大值，并用于输出。/for(k=0;k<m;k+)printf("na%d=%12.10f b%d=%12

19、.10f c%d=%12.10f",k,ak,k,bk,k,ck);getch();return 0;double fa(double x)return(1.0/5.0)*(x-2.0);/书上所给的计算要求的公式。下同。/double fb(double x)return(1.0/(3.0*x*x)*tan(x*x*x);double fc(double x)return(x-1.0)*log(x)/(log(x)+1.0-1.0/x);第四题：#include<stdio.h>#include<math.h>/方便各种数学函数的调用。/#include&l

20、t;float.h>/为了调用FLT_EPSILON，防止出现分母为零的情况。/#define N 1000/保证有足够多的运算次数能达到比较精确的结果，且便于修改。/#define pre 0.0000000001/使循环能够在达到某一误差精度后停下来。/double main()double f(double t);/声明f函数。/double r(double t);/声明r函数。/double tN,errN,relerrN,difN;/使结果能够以数组形式输出，便于分析。/char ch1="t",ch2="dif",ch3="err",ch4="relerr"/对输出进行标注。/int k,i;/定义i，以便在达到所需精度之后，对所有结果以及对最终结果进行输出。/t0=5.0,t1=15.0;/给出任意合法的初始值。/err0=fabs(t1-t0),relerr0=err0/fabs(t1+FLT_EPSILON),dif0=t1-t0;/对各个结果的初始值进行运算。/for(k=1;k<N;k+)tk+1=tk-(f(tk)*(tk-tk-1)/(f(tk)-f(tk