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

1、非线性方程的数值计算方法实验【摘要】在利用数学工具研究社会现象和自然现象，或解决工程技术等问题时，很多问题都可以归结为非线性方程的求解问题，无论在理论研究方面还是在实际应用中，求解非线性方程都占了非常重要的地位。综合当前各类非线性方程的数值解法，通过比较分析，二分法，迭代法，牛顿拉夫森方法，迭代法的收敛阶和加速收敛方法，以上的算法应用对某个具体实际问题选择相应的数值解法。关 键 词 非线性方程；二分法；迭代法；牛顿-拉夫森法；割线法等。1、 实验目的 通过本实验的学习，应掌握非线性方程的数值解法的基本思想和原理，深刻认识现实中非线性方程数值的意义；明确代数精度的概念；掌握二分法、不动点迭代法、

2、牛顿迭代法、割线法等常用的解非线性方程的方法；培养编程与上机调试能力。2、 实验原理二分法：单变量函数方程： f（x）=0其中，f(x)在闭区间a，b上连续、单调，且f(a)*f(b)<0,则有函数的介值定理可知，方程f（x）=0在（a，b）区间内有且只有一个解，二分法是通过函数在区间端点的符号来确定所在区域，将有根区间缩小到充分小，从而可以求出满足给定精度的根的近似值。下面研究二分法的几何意义： 设=1, =b, 区间，中点= 及，若=0，则，若 f()*f()<0，令=，=，则根 ,中，这样就得到长度缩小一半的有根区间,，若 f()*f()<0，令=，=，则根 ,中，这样

3、就得到长度缩小一半的有根区间,即f()f()<0,此时-=,对有根区间,重复上述步骤，即分半求中点，判断中电处符号，则可得长度有缩小一半的有根区间, 如图所示： 重复上述过程，第n步就得到根的近似序列及包含的区间套，如下：（1）（2）（3）-=(4) 且|-| (n=1,2,3.)显然lim，且以等比数列的收敛速度收敛于，因此用二分法求f（x）=0的实根可以达到任意指定精度。（二）迭代法：对给定的方程，将它转达换成等价形式：。给定初值，由此来构造迭代序列，如果迭代法收敛，即，有，则就是方程的根。在计算中当小于给定的精度控制量时，取为方程的根。（三）牛顿迭代法：设方程f(x)=0在其根的某

4、个领域U(，)内有一阶连续导数，且f() 0。求f(x)=0的根，首先要将f(x)=0转化为等价形式，并使 (x)满足不动点迭代的一般理论。 于是我们令 (x)=x+h(x)f(x),可由 ()=0来确定h(x)的结构，根据(x)=1+h()f()+h()f(x1)=1+h()f()=0可得h()=-1/f() ，由于f(x) 0，且f(x) 连续，因此当h(x)=-1/f(x) 时， h(x1)=0，即令 (x)=x-f(x)/f(x)， 从而有迭代格式 = （k=0,1,2，.）由于, , .都在U领域里，从而当B比较小时，可用f()可近似代替f(),= - ，此方法称为牛顿迭代法。（四）

5、割线法：设,为方程f(x)=0的两个近似根。用差商得：f()-f()/ - ，代替牛顿迭代公式中的导数f()，于是得到如下的迭代公式： =-。下面研究割线法的几何意义： 经过点（,f()）及点（,f()）两点作割线，其点斜式方程为： Y=f（）- ，其零点为X=- 把X用表示即得到迭代格式，需要两个初值此割线与 X 轴交点的横坐标就是新的近似值，如图所示。 三、实验内容（1） 、实验描述 1.P40.1:参照程序2.1求解出单调收敛的不动点。 2.P49.1:已知初值，时间和末值，求解汇率。 4.P69.1：已知运动方程求解运动时间和距离。 （二）、实验题目 1.使用程序2.1求解下面每个函数

6、的不动点（尽肯能多）近似值，答案精确到 小数点后12位。同时，构造每个函数的图和直线y=x来显示所有不动点。(d) 2. 如果在240个月内每月付款300美元，求解满足全部现金A为500000美元的 汇率I的近似值（精确到小数点后10位）。4.设投射方的运动程为 （a） 求当撞击地面时经过的时间，精确到小数点后10位。（b） 求水平飞行行程，精确到小数点后10位。四、实验结果及分析 1. P40.1（d）： 算法：（1）输入函数g，p0,tol,max1,令k=2。 （2）判断k>max1是否成立，如成立输出结果，如不成立，执行（3）。 （3）令p(k)=g(p(k-1),err=|p(

7、k)-p(k-1)|。 (4)判断err<tol是否成立，如成立，输出结果，如不成立，令k=k+1,跳转到（2）。startInput g,p0,tol;k=2k>max1YNLet p(k)=g(p(k-1);err=|p(k)-p(k-1)|k=k+1err<to1outputYNoutputend令，简单画出与的图像，可令x=1,。图4. 与当输入k,p,err,P=fixpt(g,1,1e-12,100)输出k =2 p =1 err =0P =1 1 通过图形我们知道不动点有2个，X=1是吸引不动点.因而函数只有一个吸引迭代点。算法：（1） 输入f,a,b,delt

8、a,max1,令k= 1。（2） f(a)=feval(f,a);f(b)=feval(f,b);k=1.（3）判断f(a)*f(b)>0是否成立，若不成立输出结果，若成立执行步骤（4）。（4）判断k>max1是否成立,若成立输出结果，若不成立执行步骤（5）。（5）令dx=yb*(b-a)/(yb-ya),c=b-x,yc=feval(f,c). (6) 判断yc=0是否成立，若成立输出结果，若不成立执行步骤（7）。（7）判断yc*yb>0是否成立，若成立，令b=c,yb=yc.若不成立,令a=c,ya=yc.( 8)令err=|b-a|.(9)判断err<delta是

9、否成立，若成立，输出结果。若不成立，令k=k+1,执行（4）。start流程图：Input f,a,b,delta,max1，k=k+1Yf(a)*f(b)>0NoutputYk>max1Ndx=yb*(b-a)/(b-a)c=b-xyc=feval(f,c)Yyc=0NYNyc*yb>0b=cyb=yca=cya=ycerr=|b-a|Nerr<deltaYoutputend当输入c,err,yc=regula(f,0.15,0.16,1e-10,100)输出err =-1.257285475730896e-008yc =2.775557561562891e-017算

10、法：（1） 输入h,r,p0,p1,delta,max1.令k=1.（2） 判断k>max1是否成立，若成立，输出；若不成立，执行（3）。（3） 令p2=p1-feval(f,p1)*(p1-p0)/(feval(f,p1)-feval(f,p0),err1=|p1-p0|,err2=|feval(r,p1)-feval(r,p0)|,p0=p1,p1=p2 (4) 判断err1<delta与err2<delta是否成立，若成立，输出；若不成立，令k=k+1,执行步骤（2）。start流程图：Input f,r,p0,p1,delta max1 k=1k>max1P2=

11、p1-feval(f,p1)*(p1-p0)/(feval(h,p1)-feval(h,p0)Err1=|p1-p0| err2=|feval(r,p1)-feval(r,p0)|C=feval(r,p1) p0=p1 p1=p2Err2<delta&err2<deltaoutputend输入p1,err1,err2,k,c=secant(h,r,9,9.5,1e-10,100)输出p1 = 9.087899668785003err2 =2.523847797419876e-011k = 6c = 1.090547960254200e+003五、实验结论1.采用的不动点迭代

12、法，只能计算出吸引不动点，不能算出排斥迭代点。 2.计算得所求的利率为。4. 计算得投射体撞击地面经过的时间p1 = 9.087899668785003飞行行程c = 1.090547960254200e+003。附件（代码）：1(d). plot the figure of x.(x-cos(x) with y=xx=-2:0.05:2.4;y=x;g=x.(x-cos(x);plot(x,g,x,y)xlabel('x');ylabel('x.(x-cos(x) y=x');%definite function g(x)function b=g(x)b=x.

13、(x-cos(x);function k,p,err,P=fixpt(g,p0,tol,max1)%input-g is iteratin function input as string 'g'% -p0 is initial guess for the fixed point% -tol is the tolerance% -max1 is the maximum number of iterations%output-k is the number of iteratins that were carried out% -p is approximation to the

14、 fixed point% -err is the error in the approximation % -P contains the sequencepnP(1)=p0;for k=2:max1 P(k)=feval(g,P(k-1); err=abs(P(k)-P(k-1); relerr=err/(abs(P(k)+eps); p=P(k); if (err<tol)|(relerr<tol),break,endendif k=max1 disp('maximum number of iterations exceeded')endP=P'2.

15、%definite the function f(x)function A=f(x)p=300;n=240;A=12*p*(1+x/12)n-1)/x-500000;function c,err,yc=regula(f,a,b,delta,max1)%input-f is the function input as a string 'f'% -a and b are the left and right points% -delta is the tolerance for the zero% -max1 is the maximum number of iterations

16、%output-c is the zero% -yc=f(c)% -err is the error estimate for cya=feval(f,a);yb=feval(f,b);if ya*yb>0 disp('ya,yb are not suitable ') returnend for k=1:max1 dx=yb*(b-a)/(yb-ya); c=b-dx; yc=feval(f,c); if yc=0,break; elseif yc*yb<0 a=c;ya=yc; else b=c;yb=yc; end err=abs(b-a); if err&l

17、t;delta break,endend4.%definite function r(t)function x=r(t)x=2400*(1-exp(-t/15);%definite function h(t)function y=h(t)y=9600*(1-exp(-t/15)-480*t;function p1,err1,err2,k,c=secant(h,r,p0,p1,delta,max1)%input-h is the function of altitude inpput as string 'h'% -r is the function of path input as a string 'r'% -p0 and p1 is the initial approximations to a zero% -delta is the tolerance for p1 and y% -max1 is the maximum number of iterations%output-p1 is the secant method approximation to the zero% -err1 is estimate for p1%