数值计算课程设计报告

数值计算方法课程设计报告课程设计名称：数值计算方法课程设计题目：非线性方程求根年级专业：组员姓名学号：指导教师：完成时间：非线性方程求根一、问题提出随着科学技术，生产力经济的发展，在科学与工程计算中存在着大量方程求根问题，例如贷款购房问题，工厂的最佳订货问题等都需要求解一类非线性方程的根，首先根据实际问题列出数学模型，确定变量，给出各个条件及相关函数；然后对建立的模型进行具体分析和研究，选择合适的求解方法；编写函数的程序，用计算机求出方程的解，通过所求解分析具体情况．求解非线性方程的问题有以下几种基本方法。二分法简单易行，但收敛较慢，仅有线性收敛速度。而且该方法不能用于求偶数重根或复根，但可以用来确定迭代法的初始值。牛顿法是方程求根中常用的一种迭代方法，它除了具有简单迭代法的优点外，还具有二阶收敛速度（在单根邻近处）的特点，但牛顿法对初始值选取比较苛刻（必须充分靠近方程的根），否则牛顿法可能不收敛。弦截法是牛顿法的一种修改，虽然比牛顿法收敛慢，但因它不需计算函数的导数，故有时宁可用弦截法而不用牛顿法，弦截法也要求初始值必须选取得充分靠近方程的根，否则也可能不收敛。二、背景分析代数方程的求根问题是一个古老的数学问题。理论上，n次代数方程在复数域内一定有n个根(考虑重数)。早在16世纪就找到了三次、四次方程的求根公式，但直到19世纪才证明大于等于5次的一般代数方程式不能用代数公式求解，而对于超越方程就复杂的多，如果有解，其解可能是一个或几个，也可能是无穷多个。一般也不存在根的解析表达式。因此需要研究数值方法求得满足一定精度要求的根的近似解。牛顿迭代法是牛顿在17世纪提出的一种求解方程．多数方程不存在求根公式，从而求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要．而在各种科学和工程计算中往往要用到非线性方程组的求解，而牛顿法又是最基础的迭代法，在各种计算力学、控制工程等领域中发挥了不可代替的作用．而在数值计算中，非线性方程组的求解同样具有重要意义．随着计算机技术的成熟和高速发展，对于非线性方程求根问题出现了大量的数学软件（如MATLAB,SAS,SPSSD等），计算机已经成为工程师应用数学解决工程问题的主要运算工具．同时，工程专业的学生对数学教育的需求重点正在从手工演绎和运算能力的培养转变到结合计算机软件进行建模、求解和论证能力的培养．我们采用Matlab数学软件平台，通过实例比较了二分法、牛顿迭代法、弦截法三种基本方法的优缺点。三、基本算法思想与实现（1）二分法单变量函数方程：f（x）=0其中，f(x)在闭区间[a，b]上连续、单调，且f(a)*f(b)0,则有函数的介值定理可知，方程f（x）=0在（a，b）区间内有且只有一个解*x，二分法是通过函数在区间端点的符号来确定*x所在区域，将有根区间缩小到充分小，从而可以求出满足给定精度的根*x的近似值。下面研究二分法的几何意义：设1a=1,1b=b,区间11,ba，中点1x=211ba及1xf，若1xf=0，则=1x，若f(1a)*f(1x)0，令2a=1a，2b=1x，则根*x[2a,2b]中，这样就得到长度缩小一半的有根区间[2a,2b]，若f(1b)*f(1x)0，令2a=1x，2b=1b，则根*x[2a,2b]中，这样就得到长度缩小一半的有根区间[2a,2b],即f(2a)f(2b)0,此时2b-2a=211ab,对有根区间[2a,2b]重复上述步骤，即分半求中点，判断中电处符号，则可得长度有缩小一半的有根区间[2a,2b],如图所示：重复上述过程，第n步就得到根*x的近似序列nx及包含*x的区间套，如下：（1）...],....[],[],[2211nnbababa（2）],[,0)()(*nnnnbaxbfaf（3）na-nb=)(1121nnba=…=12nab(4),2nnnbax且|*x-nx|12nab(n=1,2,3…..)显然lim，且以等比数列的收敛速度收敛于*x，因此用二分法求f（x）=0的实根*x可以达到任意指定精度。（2）牛顿迭代法设方程f(x)=0在其根*x的某个领域U(*x，)内有一阶连续导数，且f’(*x)≠0。求f(x)=0的根*x，首先要将f(x)=0转化为等价形式，并使(x)满足不动点迭代的一般理论。于是我们令(x)=x+h(x)f(x),可由‘(1x)=0来确定h(x)的结构，根据’(x)=1+h’(*x)f(*x)+h(*x)f’(x1)=1+h(*x)f’(*x)=0可得h(*x)=-1/f’(*x)，由于f’(x)≠0，且f’(x)连续，因此当h(x)=-1/f’(x)时，h’(x1)=0，即令(x)=x-f(x)/f‘(x)，从而有迭代格式1kx=)(')(kkkxfxfx（k=0,1,2，…..）由于1x,2x,3x…….都在U领域里，从而当B比较小时，可用f’(0x)可近似代替f’(kx),1kx=kx-)()(0xfxfk，此方法称为牛顿迭代法下面研究牛顿法的几何意义：设r是方程f（x）=0的根，选取0x作为的r初始近似值，经过(0x,f(0x))做曲线y=f(x)的切线的方程：y=f(0x)+f’(0x)(x-0x)，求出L与x的交点的横坐标1x=0x-f(0x)/f’(0x),称1x为r的一次近似值经过点（1x,f(1x)）做切线y=f(x)的切线，并求出该切线与x轴的交点横坐标：2x=1x-f(1x)/f’(1x),2x称为r的二次近似值，重复以上操作可以得到r的近似值序列。下述三个定理分别讨论了牛顿法的收敛性质：定理1：对于方程f(x)=0，设f（x）在[a，b]上有二阶连续导数且满足下述条件：（1）f(a)f(b)0;（2）f’(x)0,)(xf0,对任意的x[a,b];（3）存在0x[a,b]，使f（0x）)(0xf0，则由牛顿法产生的迭代序列nx收敛于f(x)=0的根*x，且)(2)()(**2**1limxfxfxxxxkkk定理2：对于方程f(x)=0，设f（x）在[a，b]上有二阶连续导数且满足下述条件：（1）f(a)f(b)0;（2）对任意的x[a,b]，f’(x)0,)(xf0（3）)()(afafb-a,)()(bfbfb-a则对于任何0x[a,b]，由牛顿法产生的迭代序列nx收敛于f(x)=0的根*x定理3：设*x是方程f(x)=0的根，在*x的某个开区间内)(xf连续且f’(x)0，则存在0,当0x【*x-,*x+】时，由牛顿迭代法1kx=)(')(kkkxfxfx（k=0,1,2，…..）式产生的序列nx是以不低于二阶的收敛速度收敛到*x.（3）弦截法设kx,1kx为方程f(x)=0的两个近似根。用差商得：f(kx)-f(1kx)/kx-1kx，代替牛顿迭代公式中的导数f’(kx)，于是得到如下的迭代公式：1kx=kx-)()()()(11kkkkkxxxfxfxf。下面研究割弦法的几何意义：经过点（kx,f(kx)）及点（1kx,f(1kx)）两点作割线，其点斜式方程为：Y=f（kx）-)()()(11kkkkkxxxxxfxf，其零点为X=kx-)()()()(11kkkkkxxxfxfxf把X用1kx表示即得到迭代格式，它又称为双点弦割法，需要两个初值此割线与X轴交点的横坐标就是新的近似值1kx，所以弦截法又称为割线法，如图所示。下面三个定理为弦割法收敛定理：定理1：设f（x）在其零点*x的邻域U（*x，）=[*x-，*x+](0)内有二阶连续导数，0)(*xf，则当0xU（*x，）时，由割弦法式产生的序列nx收敛于*x，且收敛的阶为。定理2：设)(xf在区间[a,b]上连续，且满足下述三点（1）f(a)f(b)0;（2）对任意的x[a,b]，有f’(x)0,)(xf0（3）)()(afafb-a,)()(bfbfb-a则对于任意初始0x，1x[a,b]，由弦割法产生的迭代序列nx收敛于f(x)=0唯一的根*x定理3：设在其零点*x的邻域U(*x,)=[*x-，*x+](0)内有二阶连续导数，f’(x)0则当0xU(*x,)时，由弦割1kx=kx-)()()()(11kkkkkxxxfxfxf式产生的序列nx收敛于*x，且收敛的阶为。四、具体应用实例分析求解033)(23xxxxf在5.1附近的根。（1）二分法建立erfen-M文件：function[k,x,wuca,yx]=erfen(a,b,abtol)a(1)=a;b(1)=b;ya=fun(a(1));yb=fun(b(1));%程序中调用的为函数ifya*yb0,disp('注意：ya*yb0,请重新调整区间端点a和b.'),returnendmax1=-1+ceil((log(b-a)-log(abtol))/log(2));%ceil是向方向取整fork=1:max1+1a;ya=fun(a);b;yb=fun(b);x=(a+b)/2;yx=fun(x);wuca=abs(b-a)/2;k=k-1;[k,a,b,x,wuca,ya,yb,yx]Ifyx==0a=x;b=x;elseifyb*yx0b=x;yb=yx;elsea=x;ya=yx;endifb-aabtol,return,endendk=max1;x;wuca;yx=fun(x);建立FUN函数文件：functiony=fun(x)y=x.^3+x.^2-3*x-3;画图：x=[-10::10];y=fun(x);plot(x,y);由图，我们选取区间[-6，6]输入程序：[k,x,wuca,yx]=erfen(-6,6,运行结果：k=13；x=；wuca=；yx=（2）牛顿迭代法建立newtonqx-M文件：function[k,xk,yk,piancha,xdpiancha]=newtonqx(x0,tol,ftol,gxmax)x(1)=x0;fori=1:gxmaxx(i+1)=x(i)-fnq(x(i))/(dfnq(x(i))+eps);piancha=abs(x(i+1)-x(i));xdpiancha=piancha/(abs(x(i+1))+eps);i=i+1;xk=x(i);yk=fnq(x(i));[(i-1)xkykpianchaxdpiancha]if(abs(yk)ftol)&((pianchatol)|(xdpianchatol))k=i-1;xk=x(i);[(i-1)xkykpianchaxdpiancha]return;endendifigxmaxdisp('请注意：迭代次数超过给定的最大值gxmax。')k=i-1;xk=x(i);[(i-1)xkykpianchaxdpiancha]return;end[(i-1),xk,yk,piancha,xdpiancha]';建立FNQ原函数文件：functiony=fnq(x)y=x.^3+x.^2-3*x-3;建立DFNQ导函数文件：functiony=dfnq(x)y=3*x.^2+2*x-3;输入程序：[k,xk,yk,piancha,xdpiancha]=newtonqx,,,20)输出结果：k=3；xk=；yk=；piancha=；xdpiancha=（3）弦截法建立gexian文件：function[k,piancha,xdpiancha,xk,yk]=gexian(x01,x02,tol,ftol,gxmax)x(1)=x01;x(2)=x02;fori=2:gxmaxu(i)=fnq(x(i))*(x(i)-x(i-1));v(i)=fnq(x(i))-fnq(x(i-1));x(i+1)=x(i)-u