下载文档原格式
西南交⼤数值分析⾮线性⽅程组的五种解法⽬录摘要 (2)1 绪论 (3)2 五种解法 (3)2.1 ⼆分法 (3)2.1.1 ⼆分法简介 (3)2.1.2⼆分法的MATLAB程序 (3)2.2 不动点迭代法（简单迭代法） (4)2.2.1 不动点迭代法简介 (4)2.2.2 不动点迭代法的MATLAB程序 (5)2.3 ⽜顿法 (5)2.3.1 ⽜顿法简介 (5)2.3.2 ⽜顿法的MATLAB程序 (5)2.4 简易⽜顿法 (6)2.4.1 简易⽜顿法简介 (6)2.4.2 简易⽜顿法的MATLAB程序 (6)2.5 割线法 (6)2.5.1 割线法简介 (6)2.5.2 割线法的MATLAB程序 (7)3 例⼦计算及⽐较分析 (7)4 结论 (11)参考⽂献 (12)摘要本论⽂介绍了⼆分法、不动点迭代法、⽜顿法、简易⽜顿法、割线法五种算法原理，然后进⾏了MATLAB编程，得到能求解⾮线性⽅程的根的程序。
本⽂分别⽤这五种⽅法的MATLAB程序对五个例⼦进⾏了计算，得到各种⽅法所需的迭代次数，迭代精度，迭代时间等，从⽽分析⽐较五种⽅法的优缺点。
关键词：⾮线性⽅程⼆分法简单迭代法⽜顿法简易⽜顿法割线法1 绪论在科学⼯作中经常出现这类问题，即求解⾮线性⽅程或⾮线性⽅程组—求x 使得f(x)=0或求X=(x1,x2,?,x n)T使得F(x)=0。
本论⽂采⽤5种⽅法即⼆分法、不动点迭代法（简单迭代法）、⽜顿法、简易⽜顿法、割线法，通过对原理的理解进⾏了MATLAB 编程，然后对⼏个例⼦进⾏各种解法计算，进⾏⽐较分析，从⽽发现各种算法的优势与不⾜，增加对各种算法的理解。
作者所使⽤的计算机配置如表1-1所⽰。
表1-1 计算平台简介2 五种解法2.1 ⼆分法2.1.1 ⼆分法简介若f是区间[a,b]的连续函数，且f (a) f (b) < 0，则f在[a,b]内必有⼀个零点。
因为f (a) f (b) < 0，所以函数f在区间[a,b]上改变符号，因此它在这个区间内⾄少存在⼀个零点。
1大学数学实验 实验报告 | 2014/4/5一、 实验目的1、学习用Matlab 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2、通过实例学习用线性代数方程组解决简化问题。
二、 实验内容项目一：种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应维持不变。
种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n ,当年年龄k 的种群数量记作x k ，繁殖率记作b k （每个雌性个体1年的繁殖的数量），自然存活率记作s k （s k =1−d k ，d k 为1年的死亡率），收获量记作ℎk ，则来年年龄k 的种群数量x ̌k 应该为x ̌k =∑b k n k=1x k , x ̌k+1=s k x k −ℎk , (k=1,2,…,n -1)。
要求各个年龄的种群数量每年维持不变就是要求使得x ̌k =x k , (k=1,2,…,n -1).(1) 如果b k , s k 已知，给定收获量ℎk ，建立求各个年龄的稳定种群数量x k 的模型（用矩阵、向量表示）.(2) 设n =5,b 1=b 2=b 5=0,b 3=5,b 4=3,s 1=s 4=0.4,s 2=s 3=0.6,如要求ℎ1~ℎ5为500,400,200,100,100，求x 1~x 5.(3) 要使ℎ1~ℎ5均为500，如何达到？问题分析：该问题属于简单的种群数量增长模型，在一定的条件（存活率，繁殖率等）下为使各年龄阶段的种群数量保持不变，各个年龄段的种群数量将会满足一定的要求，只要找到种群数量与各个参量之间的关系，建立起种群数量恒定的方程就可以求解出各年龄阶段的种群数量。
模型建立：根据题目中的信息，令x ̌k =x k ，得到方程组如下：{x ̌1=∑b k nk=1x k =x 1x ̌k+1=s k x k −ℎk =x k+1整理得到：{−x 1∑b k nk=1x k =0−x k+1+s k x k =ℎk2 大学数学实验 实验报告 | 2014/4/52写成系数矩阵的形式如下：A =[b 1−1b 2b 3s 1−100s 2−1…b n−1b n0000⋮⋱⋮000000000⋯00−10s n−1−1]令h =[0, ℎ1,ℎ2,ℎ3,…,ℎn−2,ℎn−1]Tx =[x n , x n−1,…,x 1]T则方程组化为矩阵形式：Ax =h ，即为所求模型。
实验报告一题目：非线性方程求解摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。
前言：（目的和意义）掌握二分法与Newton法的基本原理和应用。
数学原理：对于一个非线性方程的数值解法很多。
在此介绍两种最常见的方法：二分法和Newton 法。
对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。
重复运行计算，直至满足精度为止。
这就是二分法的计算思想。
Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式产生逼近解x*的迭代数列{x k}，这就是Newton法的思想。
当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。
另外，若将该迭代公式改进为其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。
程序设计：本实验采用Matlab的M文件编写。
其中待求解的方程写成function的方式，如下function y=f(x);y=-x*x-sin(x);写成如上形式即可，下面给出主程序。
二分法源程序：clear%%%给定求解区间b=1.5;a=0;%%%误差R=1;k=0;%迭代次数初值while (R>5e-6) ;c=(a+b)/2;if f12(a)*f12(c)>0;a=c;elseb=c;endR=b-a;%求出误差k=k+1;endx=c%给出解Newton法及改进的Newton法源程序：clear%%%% 输入函数f=input('请输入需要求解函数>>','s')%%%求解f(x)的导数df=diff(f);%%%改进常数或重根数miu=2;%%%初始值x0x0=input('input initial value x0>>');k=0;%迭代次数max=100;%最大迭代次数R=eval(subs(f,'x0','x'));%求解f(x0)，以确定初值x0时否就是解while (abs(R)>1e-8)x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x'));R=x1-x0;x0=x1;k=k+1;if (eval(subs(f,'x0','x'))<1e-10);breakendif k>max;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值ss=input('maybe result is error,choose a new x0,y/n?>>','s');if strcmp(ss,'y')x0=input('input initial value x0>>');k=0;elsebreakendendendk;%给出迭代次数x=x0；%给出解结果分析和讨论：1.用二分法计算方程在[1，2]内的根。
目录一、绪论------------------------------------------------------------------------------------- 2-2二、线性方程组直接解法列主元高斯LU LDL T GG T-------------------- 3-6二、线性方程组迭代法----------------------------------------------------------------- 7-10 三、四、非线性方程组数值解法二分法不动点迭代---------------------- 11-13五、非线性方程组数值解法牛顿迭代下山弦截法----------------- 14-15六、插值线性插值抛物线插值------------------------------------------------ 16-18七、插值Hermite插值分段线性插值-----------------------------------------19-22八、拟合------------------------------------------------------------------------------------ 23-24九、数值积分----------------------------------------------------------------------------- 25-29十、常微分方程数值解法梯形欧拉改进----------------------------------- 30-32 十一、常微分方程数值解法龙格库塔------------------------------------------ 33-35绪论1-1 下列各数都是经过四舍五入得到的近似值 ,试分别指出它们的绝对误差限,相对误差限和有效数字的位数.X 1 =5.420, X 2 =0.5420, X 3 =0.00542, X 4 =6000, X 5 =0.6×105注：将近似值改写为标准形式X 1 =（5*10-1+4*10-2+2*10-3+0*10-4）*101 即n=4,m=1 绝对误差限｜△X 1｜=｜X *1-X 1｜≤ 12×10m-n ＝12×10-3 相对误差限｜△r X 1｜= ｜X∗1−X1｜｜X∗1｜≤｜X∗1−X1｜｜X1｜= 12×10-3/5.4201-2 为了使101/2 的相对误差小于0.01%, 试问应取几位有效数字?1-3 求方程x 2 -56x+1=0的两个根, 使它们至少具有4位有效数字( √783≈27.982)注：原方程可改写为(x-28)2=783线性方程组解法（直接法）2-1用列主元Gauss消元法解方程组解：回代得解：X1=0 X2=-1 X3=12-2对矩阵A进行LU分解，并求解方程组Ax=b,其中解：(注：详细分解请看课本P25)A=(211132122)→(211(1/2)5/23/2(1/2)3/23/2)→(2111/25/23/21/2(3/5)3/5)即A=L×U=(11/211/23/51)×(2115/23/23/5)先用前代法解L y＝P b 其中P为单位阵（原因是A矩阵未进行行变换）即L y＝P b 等价为(11/211/23/51)(y1y2y3)＝(111)(465)解得 y 1=4 y 2=4 y 3=35再用回代解Ux ＝y ，得到结果x即Ux ＝y 等价为(2115/23/23/5)(x 1x 2x 3)＝(y 1y 2y 3)＝(443/5) 解得 x 1=1 x 2=1 x 3=1即方程组Ax=b 的解为x ＝(111)2-3 对矩阵A 进行LDL T 分解和GG T 分解，求解方程组Ax=b,其中A=(164845−48−422) ， b ＝(123)解：（注：课本 P 26 P 27 根平方法）设L=(l i j )，D=diag(d i )，对k=1,2,…,n,其中d k =a kk －∑l kj 2k−1j=1d jl ik =(a ik −∑l ij l kj k−1j=1d j )/ d k 即d 1=a 11－∑l 1j 20j=1d j =16－0=16因为 l 21=(a 21−∑l 2j l 1j 0j=1d j )/ d 1=a 21/ d 1=416=14 所以d 2=a 22－∑l 2j 21j=1d j =5－(14)2d 1=4同理可得d 3=9 即得 D=(1649)同理l 11=(a 11−∑l ij l 1j 0j=1d j )/ d 1=1616=1=l 22=l 33 l 21=(a 21−∑l 2j l 1j 0j=1d j )/ d 1=416=14 l 31=(a 31−∑l 3j l 1j 0j=1d j )/ d 1=816=12 l 32=(a 32−∑l 3j l 2j 1j=1d j )/ d 2=−4−12×14×164=−64=－32即L=(114112−321) L T=(114121−321) 即LDL T分解为A=(114112−321)(1649)(114121−321)解解：A=(164845−48−422)→(41212−32−33)故得GG T分解：A=(4122−33)(4122−33) LDL T分解为A=(114112−321)(1649)(114121−321) 由(114112−321)(y 1y 2y 3)＝(123) ，得(y 1y 2y 3)＝(0.250.8751.7083)再由(4122−33)(x 1x 2x 3)＝(0.250.8751.7083) ，得(x 1x 2x 3)＝(−0.54511.29160.5694)2-4 用追赶法求解方程组：解：(4−1−14−1−14−1−14−1−14)→(4−14−1154−415−15615−1556−120956−56209−1780209)由(4−1154−15615−120956−1780209)(y1y2y3y4y5)＝(100200)，得(y1y2y3y4y5)＝(256.66671.785700.4784753.718)再由(1−141−4151−15561−562091)(x1x2x3x4x5)＝(256.66671.785700.4784753.718)，得(x1x2x3x4x5)＝(27.0518.20525.769314.87253.718)线性方程组解法（迭代法）2-1 设线性方程组{4x 1−x 2+2x 3=1−x 1−5x 2+x 3=22x 1+x 2+6x 3=3（1） 写出Jacobi 法和SOR 法的迭代格式（分量形式） （2） 讨论这两种迭代法的收敛性（3） 取初值x (0)=(0，0，0)T ，若用Jacobi 迭代法计算时，预估误差 ｜｜x*-x (10)｜｜∞ （取三位有效数字）解：（1）Jacobi 法和SOR 法的迭代格式分别为Jacobi 法迭代格式SOR（2）因为A 是严格对角占优矩阵,但不是正定矩阵,故Jacobi 法收敛,SOR 法当0<ω≤1时收敛.⎪⎪⎪⎩⎪⎪⎪⎨⎧+--=-+-=+-=+++216131525151412141)(2)(1)1(3)(3)(1)1(2)(3)(2)1(1k k k k k k k k k x x x x x x xx x ⎪⎪⎪⎩⎪⎪⎪⎨⎧-++-=+-+-=+-+-+=++++++)216131()525151()412141()(3)1(2)1(1)(3)1(3)(3)(2)1(1)(2)1(2)(3)(2)(1)(1)1(1k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x ωωω（3）由(1)可见||B ||∞=3/4,且取x (0)=(0,0,0)T ,经计算可得x (1)=(1/4,-2/5,1/2)T ,于是||x (1)-x (0)||∞=1/2,所以有2-2 设方程组为{5x 1+2x 2+x 3=−12−x 1+4x 2+2x 3=202x 1−3x 2+10x 3=3试写出其Jacobi 分量迭代格式以及相应的迭代矩阵，并求解。
§5 非线性方程组的解法设非线性方程组为n i x x x f n i ,,2,10),,,(21 == (31)令 T n T n f f f F x x x x ),,,(,),,,(2121 == 则(31)式记为F(x)=0 (32)求x*向量，使得F(x*)≡0的问题，就是非线性方程组求解问题，称x*为方程组（32）式的解向量。
多元非线性方程组与一元非线性方程f(x)=0具有相同的形式，因此，我们可以类似地讨论它的迭代解法，如不动点迭代法、牛顿迭代法。
一、 简单迭代法（不动点迭代）：步骤：1.对于非线性方程组: F(x)=0 2.写成等价方程组:x=G(x)3.构造不动点迭代: x (k+1)=G(x (k))4.给定初始点x (0)，生成序列{x (k)} 若*)(limx x k k =∞→，则x *是方程组的解程序：function y=g(x)y(1)=0.7*sin(x(1))+0.2*cos(x(2)); y(2)=0.7*cos(x(1))-0.2*sin(x(2));%Iteratepro.mfunction y=iteratepro(x) x1=g(x); n=1;while (norm(x1-x)>=1.0e-6)&(n<=1000) x=x1; x1=g(x); n=n+1; endx1 nx0=[0.5 0.5];iteratepro(x0)二、牛顿迭代法原理：把解非线性方程的牛顿迭代法，推广到方程组牛顿迭代的核心是用线性函数近似非线性函数，逐步用线性方程的根近似非线性方程的根。
方法：对方程组F(x)=0，用Taylor 展开方法做线性近似。
设n i x f i ,,2,1),( =具有对),,2,1(n j x j =的二阶偏导数。
又设方程组F(x)=0有近似解),,,()()(2)(1)(k n k k k x x x x =将F （x ）在x (k)点作多元函数的Taylor 展开)()()()()(k k x x J x F x F -+=+高阶项 (33)其中矩阵J （Jacobi 矩阵）：为x（k ）处i f 的各个一阶偏导数矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=n n n n n n x f x f x f x f x f xf x f x f x f J212221212111牛顿迭代过程：记 )()(k k x x x-=∆ （34）并略去高阶项，得到近似方程组)()()(k k x F x J -=∆ （35）当det （J ）≠0时，解出)(k x∆，则有)1(+k x =)(k x +)(k x ∆ （36）再将)1(+k x代回（35）式，重复计算，这就是牛顿迭代过程。
牛顿迭代过程的收敛性⎩⎨⎧=+-==--=0sin 2.0cos 7.0),(0cos 2.0sin 7.0),(21y x y y x f y x x y x f解：第一次迭代： 先求出Jacobi 矩阵⎪⎪⎭⎫ ⎝⎛+-=y x y xJ cos 2.01sin 7.0sin 2.0cos 7.01 取初值5.0)0()0(==y x⎪⎪⎭⎫⎝⎛=1755.133560.0095885.038569.0Jdet(J)=0.42121≠0,而T y x F )018423.0,011114.0(),()0()0(-=，由（35）式，得方程组⎪⎪⎭⎫⎝⎛-=∆⎪⎪⎭⎫ ⎝⎛018423.001114.01755.133560.0095885.038569.0)0(x 解得：,)100145.8,0268234.0(3)0(T x-⨯=∆因此得新近似解50801.0,52682.0)1()1(==y x此时T y x F )1008.2,103.1(),(44)1()1(--⨯⨯=第二次迭代：将)1()1(,yx 再代入Jacobi 矩阵中，得方程组⎥⎥⎦⎤⎢⎢⎣⎡⨯-⨯-=∆⎪⎪⎭⎫ ⎝⎛--44)1(1008.2103.117474.135195.0097288.039491.0x 解得：T x )104689.8,100832.3(54)1(--⨯-⨯-=∆，因此迭代两次的近似解为50793.0,52651.0)1()1()2()1()1()2(=∆+==∆+=y y y x x xT y x F )104,0(),(6)2()2(-⨯=，精度已相当高。
牛顿迭代计算机算法步1，输入初始解向量ε,)0(X；步2，对i,j=1,2,…,n 计算)0(x x f ii ∂∂； 步3，求)()0(xF ；步4，求解方程组);()0()0(x F x J -=∆步5，求)0()0()1(x x x ∆+=；步6，若ε<-)0()1(max i i ix x ，则输出)1(x ，转步7否则)1()0(x x ⇐，转步2；步7，结束。
附牛顿迭代法程序： function y=fun5(x)y(1)=x(1)-0.7*sin(x(1))-0.2*cos(x(2)); y(2)=x(2)-0.7*cos(x(1))+0.2*sin(x(2)); y=[y(1) y(2)];function y=df(x)y=[1-0.7*cos(x(1)) 0.2*sin(x(2))0. 7*sin(x(1)) 1+0.2*cos(x(2))];%newtonpro.mfunction y=newtonpro(x0) x1=x0-fun5(x0)/df(x0); n=1;while(norm(x1-x0)>=1.0e-6)&(n<=1000) x0=x1;x1=x0-fun5(x0)/df(x0);n=n+1; endx1,n %输出满足精度的迭代次数与解 x0=[0.5 0.5];newtonpro(x0)三、最速下降法问题的提出：对于方程组 F(x)=0 ，构造各方程的平方和函数，即令∑==Φni i n x f x x x 1221)(),,,( （37）则求F(x)=0 的解等价于求多元函数),,,(21n x x x Φ的极小值问题。
这种无约束极值问题，可以采用最优化方法求解。
最速下降法(梯度法):1847年由Cauchy 提出，它是求解最优化问题的最古老且基本的数值方法。
后来提出的一些好的方法、都是这种方法的改进。
原理：设),,,(21n x x x Φ对各变元具有连续的偏导数，记其梯度为Tnx x x x ),,,()(21∂Φ∂∂Φ∂∂Φ∂=Φ∇ （38） 取Tn x x x x),,,()0()0(2)0(1)0( =为近似解，任取一个单位向量p 及正数t ，将)(x Φ在)0(x点进行Taylor 展开，即有)()(()()()0()0()0(t o p x t x tp x T +Φ∇+Φ=+Φ （39）其中o(t)是t 的高阶无穷小量，而)),(cos()()),(cos()()]([)0()0()0()0()0(p x x p x p x p x T Φ∇⋅Φ∇=Φ∇⋅⋅Φ∇=Φ∇ （40）1.选方向：当（p x),()0(Φ∇）=π时，即,)(/)()0()0(时x x p Φ∇Φ-∇=)()0(tp x +Φ向)()0(x Φ变化最快，就是说沿负梯度方向)(x Φ在点)0(x 附近变化最快。
最速下降法正是因此而得名。
2．选步长：从)0(x求)1(x，沿负梯度方向跨出一步：t ，即)()0()0()1(x t x x Φ∇-= （41）选步长原则：使))(()()0()0(x t xt F Φ∇-Φ=达极小，即求0t 使)(min ))(()0(0)0(t F x t x =Φ∇-Φ （42）重复上述过程，由)1(x求出新的近似解)()1(1)1()2(x t x xΦ∇-=最速下降法算法公式：k=0,1,2,…⎪⎩⎪⎨⎧Φ∇-=Φ∇-Φ=Φ∇-Φ+)44()()43())((min ))(()()()1()()()()(k k k k k k i k k k x t x xx t x x t x迭代计算，直至Φ的值降到很小，则终止计算。
【说明】最速下降法对任选的初始向量都能收敛，但收敛速度不如牛顿迭代法快。
将牛顿迭代法与最速下降法联合使用，可取长补短，效果更好。
最优化原理的使用：详见Matlab 优化工具箱Matlab 中非线性方程求根问题的实现：fsolve :Solve a system of nonlinear equations by a least squares method.【例】function F = fun7(x)F = [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))];x0=[-5; -5]; % Make a starting guess at the solution options=optimset('Display','iter');% Option to display output[x,fval] = fsolve(@fun7,x0,options) % Call optimizer %x=fsolve(@fun7,x0)x =0.56710.5671【例】fun = inline('sin(3*x)');x = fsolve(fun,[14],optimset('Display','off'));xx =1.0472 4.1888fun = inline('sin(3*x)');x = fsolve(fun,1,optimset('Display','off'));xx =1.0472【例】解方程组：x-0.7*sin(x)-0.2*cos(y)=0y-0.7*cos(x)+0.2*sin(y)=0x0=[0 0];fsolve('fun5',x0)ans =0.5265 0.5079solve:symbolic solution of algebraic equations【例】solve('sin(x) = 1')ans =1/2*pi【例】 solve('p*sin(x) = r') %chooses 'x' as the unknownans =asin(r/p)【例】方程组[x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3 = 0')x =[ 1][ 3]y =[ 1][ -3/2]【例】[x,y]=solve('x-0.7*sin(x)-0.2*cos(y)','y-0.7*cos(x)+0.2*sin(y)')x =.52652262191818418730769280519209y =.50791971903684924497183722688768。
