matlab&基本语句
matlab 基本语句
1.循环语句for
for&&i=s1:s3：s2
&&&&循环语句组
解释：首先给i赋值s1；然后，判断i是否介于s1与s2之间；如果是，则执行循环语句组,i=i+s3（否则，退出循环.）；执行完毕后，继续下一次循环。
例：求1到100的和，可以编程如下：
&&&for&i=1:1:100
&&&&&sum=sum+i
&&&这个程序也可以用while语句编程。
&&&注：for循环可以通过break语句结束整个for循环.
2.循环语句while
&&例：sum=0;i=1;
&&&&&&while(i&=100)
&&&&&&&sum=sum+i;i=i+1;
&&if(条件)
&&if(条件)
&&if(条件)
4.关系表达式：
&=,&,&,&=,&=,==(精确等于)
5.逻辑表达式：|(或),&(且)
6.[n,m]=size(A)(A为矩阵)
这样可以得到矩阵A的行和列数
n=length(A)，可以得到向量A的分量个数；如果是矩阵，则得到矩阵A的行与列数这两个数字中的最大值。
7.！后面接Dos命令可以调用运行一个dos程序。
8.常见函数：
poly():为求矩阵的特征多项式的函数，得到的为特征多项式的各个系数。如a=[1,0,0;0,2,0;0,0,3]，则poly(a)＝1&&&&-6&&&&11&&&&-6。相当于poly(a)＝1入^3+(-6)入^2+11入+(-6)。
compan():可以求矩阵的伴随矩阵.
sin()等三角函数。
MATLAB在数学建模中的应用(3)
　　一、程序设计概述
　　MATLAB所提供的程序设计语言是一种被称为第四代编程语言的高级程序设计语言，其程序简洁，可读性很强，容易调试。同时，MATLAB的编程效率比C/C++语言要高得多。
　　MATLAB编程环境有很多。常用的有：
　　1.　命令窗口
　　2.　word窗口
　　3.　M-文件编辑器，这是最好的编程环境。
　　M-文件的扩展名为“.m”。M-文件的格式分为两种：
&&①l M-脚本文件，也可称为“命令文件”。
&&②&M-函数文件。这是matlab程序设计的主流。l
保存后的文件可以随时调用。
二、MATLAB程序结构
　　按照现代程序设计的观点，任何算法功能都可以通过三种基本程序结构来实现，这三种结构是：顺序结构、选择结构和循环结构。其中顺序结构是最基本的结构，它依照语句的自然顺序逐条地执行程序的各条语句。如果要根据输入数据的实际情况进行逻辑判断，对不同的结果进行不同的处理，可以使用选择结构。如果需要反复执行某些程序段落，可以使用循环结构。
　　1&&顺序结构
顺序结构是由两个程序模块串接构成。一个程序模块是完成一项独立功能的逻辑单元，它可以是一段程序、一个函数，或者是一条语句。
　　看图可知，在顺序结构中，这两个程序模块是顺序执行的，即先执行&程序模块1&，然后执行&程序模块2&。
实现顺序结构的方法非常简单，只需将程序语句顺序排列即可。
2&&选择结构
　　在MATLAB中，选择结构可由两种语句来实现。
　　（1）&&if语句
　　if语句的最简单用法为：
　　　　if&&表达式；
　　　　　程序模块；
　　　　end
　　if语句的另一种用法为：
　　　　if&&表达式
　　　　　程序模块1
　　　　else
　　　　　程序模块2
　　　　end
　　例1&&使用if语句判断学生的成绩是否及格。
　　程序：
　　　clear&&
　　　n=input(’输入n=&’)&&
　　　m=60;
　　　if&&n＜ｍ，
r=’不及格’
　　&&&&&&
else 　　　　　　　　
r=’及格’
练习一：将例1写入M-文件编辑器，然后在command&window&调用这个程序。
当针对多个条件进行选择时，可以采用下面的格式：
　　　　if& 表达式1&
　&　程序模块1
&&elseif&&表达式2&
&&&&&&&&&&&&&&&&&程序模块2
&&&&&&&&&&
&……&&&……
&&&&&&&&&&
elseif&&表达式n
&&&&&&&&&&&&&&&&&
&&&程序模块n&
&&&&&&&&else&&
&&&&&&&&&&&程序模块n+1
　　例2&&将百分之的学生成绩转换为五分制输出。
　　程序：
　n=input(’输入n=&’)
　if&&n&=90
　chji=’优秀’
　elseif&&n&=80
　chji=’良好’
　elseif&&n&=70
　chji=’中等’
　elseif&&n&=60
　chji=’及格’
　chji=’不及格’
　　练习二：将例2写入M-文件编辑器，然后在command&window&调用这个程序。
　　（2）&&switch语句
　　switch语句可以替代多分支的if语句，而且switch语句简洁明了，可读性更好。其格式为：
　　　switch&&表达式&&&
　　　case&&数值1&&&&&
&&&&&&&&&&程序模块1
　　　case&&数值2
&&&&&&&&&&程序模块2
　　　……
　　　otherwise
　　　程序模块n
其中的otherwise模块可以省略。
　　switch语句的执行过程是：首先计算表达式的值，然后将其结果与每一个case后面的数值依次进行比较，如果相等，则执行该case的程序模块；如果都不相等，则执行otherwise模块中的语句。如图3所示。
　　例3&&用switch…case开关结构将百分制的学生成绩转换为五分制的成绩输出。
　switch&&fix(x/10)
&&&&case&&{10,9}
&&&&&&&&y=’优秀’
&&&&case&&8
&&&&&&&&y=’良好’
&&&&case&&7
&&&&&&&&y=’中等’
&&&&case&&6
&&&&&&&&y=’及格’
&&&&otherwise
&&&&&&&&y=’不及格’
　　练习三：将例3写入M-文件编辑器，然后在command&window&调用这个程序。
　　3&&循环结构
　　循环结构的流程图如图4所示它可以多次重复执行某一组语句。循环是计算机解决问题的主要手段。
　　在MATLAB中，循环结构可以由两种语句结构实现。
　　（1）&for…end 循环结构。其格式为：
&&&&&&&&&&&&&&&&&&&&&&&
for&i=V，循环体结构，end
其中V为一个行向量，循环变量i每次从V中取一个数值，执行一次循环体的内容，如此下去，直到完成V中的所有分量，就自动结束循环体的执行。
　　例4&&&&计算&&s=12+22+52。
　　程序：
　　a=[1&2&5&];&s=0;&
　　for&k=a,&
　　　　s=s+k^2;&
　　　s，&&
　　该例题只是为了说明for语句的格式，事实上，用下面的语句求和更为简单。
　　　　p=sum(a.^2)&&
　　练习四：建立MATLAB与word的连接，在新建的m-book中写入上面的程序，并用notebook菜单运行之。
　　循环结构里面还可以包含循环结构，形成多重循环。
　　例5&&&&设计一个九九乘法表。
　　程序：
　　　clear
　　　for&&i=1:9
　　　　for&&j=1:9
　　　　　a(i&,&j)=i*j;
　　　　end
&&&&&&&a，&&
练习五：①直接在命令窗编写上面的程序。
　　　　　　②试运行下面的程序，并加以分析：
　　程序：&disp(’&&&&&九九乘法表&’),for&i=1:9,
　　　　　　for&j=1:i,p{i}(j)=j*i;end,disp(p{i}),end&&
　　（2）&while…end
循环结构。其格式为：
　　　　　　　　while&&（表达式）
　　　　　　　　　　循环结构体
　　　　　　　　end
　　例6&&求MATLAB的相对精度
　　解：解题的思路是，让y值不断减小，直到MATLAB分不出1+y与1的差别为止。
　　程序：
　　　　y=1;&while&&1+y&1;&&y1=y;&&y=y/2;&end,y1
　　说明：
　　①for循环与while循环的区别是，for语句的循环次数是确定的，而while语句的循环次数是不确定的。
　　②一定要注意在循环结构体内设置“修改条件表达式的语句”，以免进入“死循环”。
　　③一旦出现死循环，在命令窗用ctrl+c可使程序中止。
　　④注意程序的可读性。
　　练习六：在M-文件编辑器内编写上面的脚本文件，并利用菜单或快捷按钮运行之。
　　三、M-文件
　　M文件是包含MATLAB代码的文件。M文件按其内容和功能可分为脚本M文件和函数M文件这两大类。
　　1、脚本M文件
　　脚本M文件是许多MATLAB代码按顺序组成的命令集合，不接受参数的输入和输出，与MATLAB工作区共享变量空间。脚本文件一般用来实现一个相对独立的功能，比如对某个数据集进行某种分析、绘图，求解方程等等。
前面的几个例题都是脚本文件的类型。
　　2、函数M文件
　　MATLAB的M-函数是由function语句引导的，其基本格式如下：
　　　　　　　function&[输出形参列表]&=&函数名&(输入形参列表)
　　　　　　　注释说明语句段，由%引导;
　　　　　　　函数体语句
　　函数文件需要在M-文件编辑器中编写。写完以后，按照当前的搜索路径保存。以后就可以随时调用这个函数了。
　　与脚本M文件不同的是，调用函数M文件时需要输入自变量的实际值。
随便打开一个M-文件看一看：
　　　　　　　　open&lsqlin
　　编程实例
　　例7.&定义一个函数&&f(x)=[sin(x)]2,&其中x以“度”为单位。然后再调用该函数。
　　解：在M-文件编辑器内写入下面的文件：
　　　　function&y=sinsd(x)
　　　　%自变量x以“度”为单位。
　　　　%白城师院
　　　　%数学建模协会，
　　　　a=x/180*
　　　　y=sin(a).^2;
　　将上面的文件保存之后即可随时调用。
　　　　t=sinsd(45)&
&&&　　&0.5000&&
　　　　help&sinsd&&
　　　　t=sinsd([45,60])&&
&　　&&&0.5000&&&&0.7500&&
表明该函数对元素群运算有效。上面的语句相当于
　　　　x=[45,60];t=sinsd(x)&&
&&&&　　0.5000&&&&0.7500&&
　　　　x=[45,60;30,90];f=sinsd(x)&&
&　　&&&0.5000&&&&0.7500
　　&&&&0.2500&&&&1.0000&&
　　将源文件中最后一行的“.”去掉，再运行以上两条命令，我们会发现什么？&
　　例8.&在MATLAB中，一个函数可以调用其它函数，也可以调用自身，即递归调用。下面利用递归算法编写一个函数，用来计算Fibonacci数列的第k项。
　　Fibonacci数列：
　　　　1，1，2，3，5，8，13，21，……
　　M-函数文件：
　　　　function&a=my_fibo(k)
　　　　if&k==1|k==2,a=1;
　　　　else,a=my_fibo(k-1)+my_fibo(k-2);end&&
　　将这个文件写入M-文件编辑器并以名称“my_fibo.m”保存，以后就可以调用这个函数。
　　递归算法无疑是解决某一类问题的有效方法，但不宜滥用，因为它的运算速度往往很慢。
　　　　tic,&&n=&my_fibo(26)&&,toc&&
&&&&&　　121393
　　elapsed_time&=
&　　　&&34.4290&&
　　下面我们尝试用一般的循环语句来求解这个问题。
　　　　tic,n=[1,1];for&k=3:100,n(k)=n(k-1)+n(k-2);end,toc,&&
　　elapsed_time&=
&　　　　　　&&&0.0100&&
　　　　n(1:26)&&
&&Columns&1&through&8&&
&&&&&&&1&&&&&&&1&&&&&&&2&&&&&&&3&&&&&&&5　　　&8　　　&13&&&&&&&&21
&&Columns&9&through&16&&
&&&&&&34&&&&&&55&&&&&&89&&&&&144&&&&&233　　&377&&&&&&610&&&&&&&987
&&Columns&17&through&24&&
&&&&1597&&&&2584&&&&4181&&&&6765&&&10946　
17711&&&&28657&&&&&46368
&&Columns&25&through&26&&
&&&75025&&121393&&
　　例9.&可变输入变量个数的函数
　　MATLAB提供的conv(&)函数可用来求两个多项式的乘积。对于多个多项式的连乘，则不能直接使用此函数，需要用该函数嵌套使用，用起来很不方便。下面编写一个MATLAB函数，使它能直接处理任意多个多项式的乘积问题。
　　M-函数文件:
　　　　function&a=convs(varargin)
　　　　　　a=1;
　　　　for&i=1:length(varargin)
&&&　　　　&a=conv(a,varargin{i});
　　　　end
　　形参varargin是一个特殊的字符串，它把输入变量列表转换成一个元胞数组，每一个输入变量都是这个元胞数组的一个元素。下面调用这个函数，求解
　　　　　　　　　d=(x4+2x3+4x2+5)(x+2)(x2+2x+3)
　　　　p=[1,2,4,0,5];q=[1,2];f=[1,2,3];
　　　　d=convs(p,q,f)&&
&&&&　&1&&&&&6&&&&19&&&&36&&&&45&&&&44&&&&35&&&&30&&
　　　　convs(p,q,f,[1,1],[1,3],[1,1])&&
&&&　　&1&&&&11&&&&56&&&176&&&376&&&578&&&678&&&648&&&527&&&315&90&&
　　例10&关于break、continue、return的用法。
　　当程序运行过程中出现return命令时，程序停止运行。break、continue用在循环语句中。在循环语句中，如果遇到break命令时，程序结束当前的“for”或“while”循环，转而执行它下面最近的end以下的语句；遇到continue时，跳过当次循环而继续下一次的循环，例如，原定要循环5次，但在进入第3次循环时遇到了continue，则第3次的循环被跳过，而继续第4次、第5次的循环。循环体实际上只重复执行了4次。
&&& 程序：
　　　　clear
　　　　str=’MATLAB&R14.3&version’;
　　　　for&i=1:length(str)
&&&&&if&(~isletter(str(i)))
　　　&&&&&
　&&&&continue
&& result(i)=str(i);
&result&=&
&&&&&&&&&&&&&
MATLAB&R&&&&&version&&&&&
运行这个程序并观察结果。然后再将程序中的语句continue改为break或return，看运行结果有何变化。
已投稿到：
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。《用matlab求解非线性方程组的几种方法之程序》 www.wenku1.com
用matlab求解非线性方程组的几种方法之程序日期：
第二章 非线性方程（组）的数值解法的MATLAB程序本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（Müller）法和迭代法的加速等）及其MATLAB程序，求解非线性方程组的方法及其MATLAB程序. 2.1 方程2.1 方程（方程（组）的根及其MATLAB命令 命令 2.1.2 求解方程2.1.2 求解方程（求解方程（组）的solve命令 求方程f(x)=q(x)的根可以用MATLAB命令:>> x=solve('方程f(x)=q(x)',’待求符号变量x’)求方程组fi(x1,…,xn)=qi(x1,…,xn) (i=1,2,…,n)的根可以用MATLAB命令:>>E1=sym('方程f1(x1,…,xn)=q1(x1,…,xn)');…………………………………………………….En=sym('方程fn(x1,…,xn)=qn(x1,…,xn)');[x1,x2,…,xn]=solve(E1,E2,…,En, x1,…,xn) 2.1.3 求解2.1.3 求解多项式求解多项式方程多项式方程（方程（组）的roots命令 如果f(x)为多项式，则可分别用如下命令求方程f(x)=0的根，或求导数f'(x)（见表 2-1）.表 2-1 求解多项式方程（组）的roots命令命 令xk =roots(fa) 功 能 输入多项式f(x)的系数fa（按降幂排列），运行后输出xk为f(x)=0的全部根.输入多项式f(x)的系数fa（按降幂排列），运行后输出dfa为多项式f(x)的导数f(x)的系数.输入多项式f(x)的导数f(x)的系数dfa（按降幂排列），运行后输出dfx为多项式f(x)的导数f(x).'''dfa=polyder(fa) dfx=poly2sym(dfa)2.1.4 2.1.4 求解方程（求解方程（组）的fsolve命令 如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots命令.如果非线性方程（组）是含有超越函数，则无法使用roots命令，需要调用MATLAB系统中提供的另一个程序fsolve来求解.当然，程序fsolve也可以用于多项式方程（组），但是它的计算量明显比roots命令的大.fsolve命令使用最小二乘法（least squares method）解非线性方程（组）F(X)=0的数值解，其中X和F（X）可以是向量或矩阵.此种方法需要尝试着输入解X的初始值（向量或矩阵）X0，即使程序中的迭代序列收敛，也不一定收敛到F(X)=0的根（见例2.1.8）.fsolve的调用格式： X=fsolve(F,X0)0F(X)=0，运行后输出F(X)=0解的估计值（向量或矩阵）X.要了解更多的调用格式和功能请输入:help fsolve,查看说明. 2.2 2.2 搜索根的方法及其MATLAB程序 程序 求解非线性方程根的近似值时，首先需要判断方程有没有根？如果有根，有几个根？如果有根，需要搜索根所在的区间或确定根的初始近似值（简称初始值）.搜索根的近似位置的常用方法有三种:作图法、逐步搜索法和二分法等，使用这些方法的前提是高等数学中的零点定理. 2.2.1 2.2.1 作图法及其 作图法及其MATLAB程序 程序 作函数的图形的方法很多，如用计算机软件的图形功能画图，或用高等数学中应用导数作图，或用初等数学的函数叠加法作图等.下面介绍两种作图程序.作函数y=f(x)在区间 在区间 [a,b] 的图形的] 的图形的MATLAB程序一 程序一x=a:h:b; % h是步长y=f(x); plot(x,y)grid, gtext('y=f(x)')说明：⑴ 此程序在MATLAB的工作区输入，运行后即可出现函数y=f(x)的图形.此图形与x轴交点的横坐标即为所要求的根的近似值.⑵ 区间[a,b] 的两个端点的距离 b-a 和步长h的绝对值越小，图形越精确.作函数y=f(x)在区间 在区间 [a,b]上的图形的MATLAB程序二 程序二将y=f(x)化为h(x)=g(x)，其中h(x)和g(x)是两个相等的简单函数x=a:h:b; y1=h(x); y2=g(x);plot(x, y1, x, y2)grid,gtext(' y1=h(x),y2=g(x)')说明：此程序在MATLAB的工作区输入，运行后即可出现函数y1=h(x)和y2=g(x)的图形.两图形交点的横坐标即为所要求的根的近似值.下面举例说明如何用计算机软件MATLAB的图形功能作图. 2.2.2 2.2.2 逐步搜索法及其逐步搜索法及其MATLAB程序 程序 逐步搜索法也称试算法.它是求方程f(x)=0根的近似值位置的一种常用的方法. 逐步搜索法依赖于寻找连续函数f(x)满足f(a)与f(b)异号的区间[a,b].一旦找到区间，无论区间多大，通过某种方法总会找到一个根.MATLAB的库函数中没有逐步搜索法的程序，现提供根据逐步搜索法的计算步骤和它的收敛判定准则编写其主程序，命名为zhubuss.m.逐步搜索法的MATLAB主程序 主程序输入区间端点a和b的值，步长h和精度tol，运行后输出迭代次数 k=(b-a)/h+1,方程 f(x)=0根的近似值r.function [k,r]=zhubuss(a,b,h,tol)% 输入的量--- a和b是闭区间[a,b]的左、右端点；%---h是步长；%---tol是预先给定的精度.% 运行后输出的量---k是搜索点的个数；% --- r是方程 在[a,b]上的实根的近似值，其精度是tol；X=a:h:b;Y=funs(X);n=(b-a)/h+1;m=0;X(n+1)=X(n);Y(n+1)=Y(n);for X(k)=a+k*h;Y(k)=funs(X(k)); %程序中调用的funs.m为函数sk=Y(k)*Y(k-1);if skm=m+1;r(m)=X(k);endxielv=(Y(k+1)-Y(k))*(Y(k)-Y(k-1));if (abs(Y(k))m=m+1;r(m)=X(k);endend 例2.2.4 用逐步搜索法的MATLAB程序分别求方程2x+2x-3x-3=0和sin(cos2x3)=0在区间[-2,2]上的根的近似值，要求精度是0.000 1.解 在MATLAB工作窗口输入如下程序>> [k,r]=zhubuss(-2,2,0.001,0.0001)运行后输出的结果k =4001r = -1.0 -1.0 1.2250即搜索点的个数为k =4 001，其中有5个是方程⑴的近似根，即r = -1.224 0，-1.000 0，-1.000 0，-0.999 0，1.225 0，其精度为0.000 1.在程序中将y=2.*x.^3+2.*x.^2-3.*x-3用y=sin(cos(2.*x.^3)) 代替，可得到方程⑵在区间[-2,2]上的根的近似值如下r = -1.0 -1.0 -0.0 1.33101.0 1.9200如果读者分别将方程⑴的结果与例2.2.3比较，方程⑵的结果与例2.1.2比较，将会发现逐步搜索法的MATLAB程序的优点.如果精度要求比较高，用这种逐步搜索法是不合算的.322.3 2.3 二分法及其MATLAB程序 程序 2.3.1 2.3.1 二分法二分法也称逐次分半法.它的基本思想是：先确定方程f(x)=0含根的区间（a,b），再把区间逐次二等分.我们可以根据式（2.3b）、区间[a,b]和误差ε，编写二分法求方程根的迭代次数的MATLAB命令.已知闭区间[a,b]和误差ε.用二分法求方程误差不大于ε的根的迭代次数k的MATLAB命令k=-1+ceil((log(b-a)- log(abtol))/ log(2)) % ceil是向+∞方向取整，abtol是误差ε. 2.3.2 2.3.2 二分法的MATLAB程序 程序 *二分法需自行编制程序，现提供用二分法求方程f(x)=0的根x的近似值xk的步骤和式（2.3a）编写一个名为erfen.m的二分法的MATLAB主程序如下.二分法的MATLAB主程序 主程序求解方程f(x)=0在开区间（a,b）内的一个根的前提条件是f(x)在闭区间[a,b]上连续，且f(a)?f(b)输入的量：a和b是闭区间[a,b]的左、右端点，abtol是预先给定的精度.运行后输出的量：k是使用二分法的次数.x是方程 在（a,b）内的实根x*的近似值，kk的近似值x的绝对精度限,满足wuca≤ abtol.yx=f(xk) ，即方程f(x)=0在实根x*的近似值x处的函数值.function [k,x,wuca,yx]=erfen(a,b,abtol)a(1)=a; b(1)=b;ya=fun(a(1)); yb=fun(b(1)); %程序中调用的fun.m 为函数if ya* yb>0,disp('注意：ya*yb>0,请重新调整区间端点a和b.'), returnendmax1=-1+ceil((log(b-a)- log(abtol))/ log(2)); % ceil是向+∞ 方向取整for k=1: max1+1a;ya=fun(a);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]if yx==0a=x; b=x;elseif yb*yx>0b=x;yb=elsea=x; ya=endif b-aendk=max1; yx=fun(x); 3例2.3.3 2.3.3 确定方程x-x+4=0的实根的分布情况，并用二分法求在开区间 (-2,-1)内的实根的近似值，要求精度为0.001. 解 （1）先用两种方法确定方程x3-x+4=0的实根的分布情况。的实根的分布情况。方法1 作图法.作图法.在MATLAB工作窗口输入如下程序>>x=-4:0.1:4;y=x.^3-x +4; plot(x,y)grid,gtext('y=x^3-x+4')画出函数f(x)=x3-x+4的图像，如图2-7.从图像可以看出，此曲线有两个驻点±3都在x3轴的上方，在(-2,-1)内曲线与x轴只有一个交点，则该方程有唯一一个实根，且在 (-2,-1)内.方法2 试算法 试算法.试算法.图2-7 在MATLAB工作窗口输入程序>>x=-4: 1:4,y=x.^3-x +4运行后输出结果x = -4 -3 -2 -1 0 1 2 3 4 y = -56 -20 -2 4 4 4 10 28 64由于连续函数f(x)满足f(-2)?f(-1)（2） 用二分法的主程序计算 用二分法的主程序计算.用二分法的主程序计算.在MATLAB工作窗口输入程序>> [k,x,wuca,yx]=erfen (-2,-1,0.001)运行后屏幕显示用二分法计算过程被列入表 2-3，其余结果为wuca = 9.，yx = 0.00372.4 迭代法及其2.4 迭代法及其MATLAB程序 程序 确定根的近似位置以后，接下来的工作就是将根精确化，即按某种方法将初始近似值逐步精确化，直到满足所要求的精确度为止. 上节介绍的二分法是将根精确化的方法之一，但是它的收敛速度较慢，且不能求出偶重根.迭代法可以克服这种缺陷.迭代法是求解方程的根、线性和非线性方程组的解的基本而重要的方法. 2.4.2 迭代法的2.4.2 迭代法的MATLAB程序1 迭代法需要自行编制程序.下面提供的迭代法的MATLAB程序1使用时只需输入迭代初始值x0、迭代次数k、迭代公式xk+1=?(xk)和一条命令，运行后就可以输出求迭代序列{xk}、迭代k次得到的迭代值xk和相邻两次迭代的偏差piancha =|xk-xk-1| （简称偏差）和偏差的相对误差xdpiancha=|xk-xk-1|xk的值，并且具有警报功能（若迭代序列发散，则提示用户“请重新输入新的迭代公式”；若迭代序列收敛，则屏幕会出现“祝贺您！此迭代序列收敛，且收敛速度较快”）.我们可以用这个程序来判断迭代序列的敛散性，也可以用于比较由一个方程得到的几个迭代公式的敛散性的优劣.function [k,piancha,xdpiancha,xk]=diedai1(x0,k)% 输入的量--x0是初始值,k是迭代次数x(1)=x0;for i=1:kx(i+1)=fun1(x(i));%程序中调用的fun1.m为函数y=φ(x)piancha= abs(x(i+1)-x(i)); xdpiancha= piancha/( abs(x(i+1))+eps);i=i+1;xk=x(i);[(i-1) piancha xdpiancha xk]endif (piancha >1)&(xdpiancha>0.5)&(k>3)disp('请用户注意：此迭代序列发散，请重新输入新的迭代公式')endif (piancha 3)disp('祝贺您！此迭代序列收敛，且收敛速度较快')endp=[(i-1) piancha xdpiancha xk]'; 例2.4.1 求方程f(x)=x2+2x-10的一个正根.解 在MATLAB工作窗口输入程序>> [k,piancha,xdpiancha,xk]= diedai1(2,5)运行后输出用迭代公式xk+1=(10-xk)/2的结果21.00 1.00 0.33 3.002.00 2.00 5.00 0.003.00 4.00 0.90 4.004.00 11.00 1.51 -6.005.00 11.13 0.97 -18.13请用户注意：此迭代序列发散，请重新输入新的迭代公式k=5，piancha = 11.13，xdpiancha = 0.97，xk = -18.13*由以上运行后输出的迭代序列与根x=2.316 624 790 355 40相差越来越大，即迭代序列{xk}发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于0.5.用迭代公式xk+1=10/(xk+2)运行后输出的结果[k,piancha,xdpiancha,xk]=1.00 0.00 0.00 2.002.00 0.78 0.00 2.223.00 0.36 0.73 2.584.00 0.55 0.16 2.025.00 0.28 0.15 2.30k =5,piancha =0.28,xdpiancha = 0.15,xk = 2.30可见，偏差piancha和偏差的偏差的相对误差xdpiancha的值逐渐变小，且第5次的迭代值xk =2.331 460 674 157 30与根x=2.316 624 790 355 40接近，则迭代序列{xk}收敛，但收敛速度较慢，此迭代法较为成功.用迭代公式xk+1=xk-(xk+2xk-10)/(2xk+2)运行后输出的结果[k,piancha,xdpiancha,xk]= 1.00 0.33 0.14 2.332.00 0.67 0.32 2.67 3.00 0.90 0.22 2.77 4.00 0.37 0.12 2.40 5.00 0 0 2.40祝贺您！此迭代序列收敛，且收敛速度较快k = 5; piancha = 0; xdpiancha = 0; y = 2.40. *2可见，偏差piancha和偏差的相对误差xdpiancha的值越来越小，且第5次的迭代值xk=2.316 624 790 355 40与根x=2.316 624 790 355 40相差无几，则迭代序列{xk}收敛，且收敛速度很快，此迭代法成功. * 2.4.5 迭代法的2.4.5 迭代法的MATLAB程序2function [k,piancha,xdpiancha,xk,yk]=diedai2(x0,tol,ddmax)x(1)=x0;for i=1: ddmaxx(i+1)=fun(x(i));piancha=abs(x(i+1)-x(i));xdpiancha=piancha/( abs(x(i+1))+eps);i=i+1;xk=x(i);yk=fun(x(i)); [(i-1) piancha xdpiancha xk yk]if (pianchak=i-1; xk=x(i);returnendendif i>ddmaxdisp('迭代次数超过给定的最大值ddmax')endP=[(i-1),piancha,xdpiancha,xk,yk]'; 例2.4.5 求x5-3x+1=0在0.3附近的根，精确到4位小数.解 ⑴ 构造迭代公式xk=?(xk).改写原方程为等价方程x=(x式5xk+1=(xk+1)/3 (k=0,1,2,L). 5+1)/3.这时?(x)=(x5+1)/3，初始值为x0=0.5, 迭代公⑵利用迭代法的利用迭代法的MATLAB迭代法的MATLAB程序MATLAB程序2程序2计算精确到4计算精确到4位小数的根的近似值位小数的根的近似值 y 和迭代次数k.在MATLAB工作窗口输入程序>> [k,piancha,xdpiancha,xk,yk]=diedai2(0.3,1e-4,100)运行后输出的结果[k,piancha,xdpiancha,xk,yk] =0 0.18 0.141 0.18 0.722 0.73 0.733 0.04 0.73k =3;piancha =1.697e-005;xdpiancha =3.680e-005;xk =0.3347; yk =0.3347. 2.5 迭代过程的加2.5 迭代过程的加速方法及其迭代过程的加速方法及其MATLAB程序 程序 对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度，但有时迭代过程收敛缓慢，从而使计算量变得很大，因此，迭代过程的加速是个重要的问题. 2.5.2 加权迭代法的2.5.2 加权迭代法的MATLAB程序 程序 加权迭代法的MATLAB主程序 主程序已知初始值x0、精度tol和迭代公式xk+1=?(xk)，求满足精度tol的根的近似根xk和它的函数值yk及迭代次数k.根据加权迭代的公式（2.10）和已知条件，现提供名为jasudd.m的M文件如下：function [k,xk,yk]=jasudd (x0,tol,L,ddmax)x1(1)=x0;for i=1: ddmaxx(i+1)=fun(x1(i));x1(i+1)= x(i+1)+L*( x(i+1)- x1(i))/(1-L);piancha=abs(x1(i+1)-x1(i));xdpiancha= piancha/( abs(x1(i+1))+eps);i=i+1;xk=x1(i);yk=fun(x1(i));[(i-1) xk yk]if (pianchak=i-1; xk=x1(i);endendif i>ddmax'迭代次数超过给定的最大值ddmax'endP=[(i-1),xk,yk]';例2.5.1 用（2.10）式求2e-x=0在0.85附近的一个近似根，要求精度ε=10解 在MATLAB工作窗口输入程序>> [k,xk,yk]=jasudd (0.85,1e-6,-0.855,100)运行后输出结果[k,xk,yk] =1.00 0.41 0.612.00 0.91 0.363.00 0.99 0.55k =3； xk =0.852606； yk = 0.852606. -x-6. 2.5.4 艾特肯2.5.4 艾特肯(Aitken)艾特肯(Aitken)加速方法的(Aitken)加速方法的MATLAB程序 程序 艾特肯加速方法的MATLAB程序 程序已知初始值x0、精度tol和迭代公式xk+1=?(xk)，求满足精度tol的根的近似根xk和它的函数值yk及迭代次数k，p=[k',x1',x2',x'].根据艾特肯加速方法的公式（2.11）和已知条件，现提供名为Aitken.m的M文件如下：function [k,xk,yk,p]= Aitken (x0,tol, ddmax)x(1)=x0;for i=1: ddmaxx1(i+1)=fun(x(i)); x2(i+1)=fun(x1(i+1));x(i+1)=x2(i+1)-(x2(i+1)-x1(i+1))^2/(x2(i+1)-2*x1(i+1)+ x(i));piancha=abs(x(i+1)-x(i));xdpiancha= piancha/( abs(x(i+1))+eps);i=i+1; xk=x(i);yk=fun(x(i));if (pianchak=i-1; xk=x1(i);endendif i>ddmaxdisp('迭代次数超过给定的最大值ddmax')endm=[0,1:i-1]; p=[m',x1',x2',x']; -x例2.5.3 用艾特肯加速方法求2e-x=0在0.85附近的一个近似根，要求精度ε=10-6.解 在MATLAB工作窗口输入程序>> [k,xk,yk,p]= Aitken(0.85,1e-6, 100)运行后输出结果k=3，xk=0.43，yk=0.73p = 0 0 0 0.00 1.00 0.45 0.84 0.07 2.00 0.11 0.26 0.07 3.00 0.43 0.98 0.73 2.6 牛顿2.6 牛顿(Newton)牛顿(Newton)切线(Newton)切线法及其切线法及其MATLAB程序 程序 的供程序：function [y,f]=newjushou(x)f=fnq(x); fz=fnq(x)*ddfnq(x)/((dfnq(x))^2+eps);y=abs(fz);if (ydisp('恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值f如下')elsedisp('请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值')endP=[y,f]'; 例 2.6.2 用牛顿切线法的局部收敛性判别方程 exsinx=4的近似根时，由下列初始值x0产生的迭代序列是否收敛？⑴x0=-1; ⑵x0=0; ⑶x0=1; ⑷x0=2; ⑸ x0=5.5；⑹x0=8.解 在MATLAB工作窗口输入程序>> [y,f]=newjushou(-1)运行后输出结果请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值y =139.5644f =4.3096说明：实际上，这个初始值所产生的迭代序列{xk}发散,见例 2.6.6.（2）输入程序>> [y,f]=newjushou(0)运行后输出结果说明：实际上，这个初始值所产生的迭代序列{xk}收敛到方程的根2.925 32，而不是收敛到离它最近的根1.400 81,即不是局部收敛的，见例 2.6.6.（3）输入程序>> [y,f]=newjushou(1)请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值 y = 8.0000 f = 4恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值f如下y =0.3566f =1.712681,即{xk}说明：实际上，这个初始值所产生的迭代序列{xk}收敛到离它最近的根1.400 说明：是局部收敛的，见例 2.6.6.（4）输入程序：>> [y,f]=newjushou(2)运行后输出结果请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值y =1.2593f =-2.7188说明：说明：虽然y =1.259 3>1，不满足牛顿切线法局部收敛的条件′(x)}是局部收际上，这个初始值所产生的迭代序列{x}收敛到离它最近的根1.400 81,即{xkk*敛的.这说明?′(x0)（5）输入程序>> [y,f]=newjushou(5.5)运行后输出结果请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值y =1.f =176.6400说明：实际上，这个初始值所产生的迭代序列{xk}收敛到235.619,但不是此方程的根，说明：比较初始值5.5与迭代值235.619，可见{xk}不是局部收敛的，见例 2.6.6.（6）输入程序>> [y,f]=newjushou(8)运行后输出结果恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值f如下y =0.4038f =-2.说明：说明：实际上，这个初始值所产生的迭代序列{xk}收敛到方程的根6.290 60，而不是收敛到离它最近的根9.424 46,即不是局部收敛的，见例 2.6.6. 2.6.3 牛顿切线2.6.3 牛顿切线法的牛顿切线法的MATLAB程序 程序 牛顿切线法求方程f(x)=0的近似根function [k,xk,yk,piancha,xdpiancha]=newtonqx(x0,tol,ftol,gxmax)x(1)=x0;for i=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) xk yk piancha xdpiancha] if (abs(yk)if i>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax。') end[(i-1),xk,yk,piancha,xdpiancha]'; 例 2.6.3 用牛顿切线法求方程2x3-3x2+1=0在x0=-0.4和0.9的近似根，要求精度ε=10.然后判断此方法分别在方程的根x=-0.5和1处的收敛速度.解 （1）先求方程的近似根.先求方程的近似根. 在MATLAB工作窗口输入程序>> [k,xk,yk,piancha,xdpiancha]=newtonqx(-0.4,0.001, 0.001,100)运行后输出初始值x0=-0.4的迭代结果如表 2-10所示.-3*kk =1.712e-007和相对误差xdpiancha =3.423e-007.如果将步骤④命令中的-0.4改作0.9，则运行后输出初始值x0=-0.4的迭代结果为 k =7；piancha =7.290e-004；xdpiancha =7.295e-004；xk = 0.999；yk =1.590e-006.（2）再讨论收敛速度比较初始值分别是x0=-0.4和0.9的两组结果，在x0=-0.4处迭代3次就得到单根x*=-0.5的近似值xk =-0.500；而在x0=0.9处迭代7次才得到二重根x*=1的近似值xk=0.999.可见，牛顿切线迭代法在单根处的迭代速度比二重根处的迭代速度快很多.这正如前面讨论的结果一样，即若x是f(x)=0的单根，则牛顿切线法是平方收敛的；若x是则牛顿切线法是线性收敛的；我们也可以用阶的定义判断它们的收敛速f(x)=0的二重根，度.留给读者证明.** 2.6.4 求2.6.4 求c的方法及其MATLAB程序 程序 用c的n次根function [k,xk,yk,piancha,xdpiancha,P]=kainfang(x0,c,n,tol, gxmax)x(1)=x0; for i=1: gxmaxu(i)= (x(i)^n-c)/(n*x(i)^(n-1)); x(i+1)= x(i)-u(i); 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),xk,yk,piancha,xdpiancha]if (pianchaend endif i>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax.')[(i-1),xk,yk,piancha,xdpiancha] endP=[(i-1),xk,yk,piancha,xdpiancha]';例2.6.5 求，要求精度为10解 本题介绍四种解法.-6.方法1 用求1 用求c的n次根c（当n是偶数时，是偶数时，c>0）的MATLAB程序计算.程序计算.取初始值x0=10，根指数n=2，被开方数c=113，近似根的精度tol=10-5，迭代的最大次数gxmax=100.在工作区间输入程序>> [k,xk,yk,piancha,xdpiancha]=kainfang(10,113,2,1e-5,100)运行后输出结果k=4,piancha=1.147e-011, xdpiancha=1.706e-012xk =10.65,yk =1.509e+009 可见,≈10.630 15，满足精度10.方法2 用牛顿迭代公式2 用牛顿迭代公式（用牛顿迭代公式（2.12）2.12）计算.计算.2则x-113=0，记f(x)=x2-113，f'(x)=2x.由牛顿迭代公式得，设x=113，-52xk-1131113xk+1=xk-，即xk+1=xk+(其中k=0,1，2，3，L)2xk2xk取初始值x0=10，计算结果列入表 2-12.因为，迭代次数k=4时，偏差-8-5x4-x3≈10.630 15.方法3 用牛顿切线法的MATLAB主程序计算.分别建立名为fnq.m和dfnq.m的M文件function y=fnq(x)function y=dfnq(x)y=2*x;在MATLAB工作窗口输入程序>> [k,xk,yk,piancha,xdpiancha]=newtonqx(10, 1e-5, 1e-5,100) 运行后，将输出的结果列入下表 2-13. 迭代k=4次，得到精度为10的结果≈ 10.630 15.表 2-13-5k 1 2 3 4 piancha 0.650 000 0.019 836 0.000 019 0.000 000 xdpiancha 0.061 033 0.001 866 0.000 002 0.000 000 xk10.650 000 10.630 164 10.630 146 10.630 146 yk0.422 500 0.000 393 0.000 000 0.000 000方法4 在MATLAB工作空间输入程序>> 113^(1/2)运行后输出ans = 10.65 经过四舍五入后，得到精度为10-5的结果≈10.630 15. 2.6.6 牛顿切线2.6.6 牛顿切线法的加速及其两种牛顿切线法的加速及其两种MATLAB程序 程序 当x是方程f(x)=0的m(m≥2)重根时，由牛顿切线公式（2.12）产生的迭代序列{xk}收敛到x的速度较慢，是线性收敛的.我们要提高收敛速度，可以有选择地采用已知根的重数和未知根的重数两种求重根的修正牛顿切线公式.**function [k,piancha,xdpiancha,xk,yk]=newtonxz(m,x0,tol,ftol,gxmax) x(1)=x0;for i=1: gxmaxx(i+1)=x(i)-m*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));if ((pianchak=i-1; xk=x(i);[(i-1) piancha xdpiancha xk yk] end endif i>gxmax'请注意：迭代次数超过给定的最大值gxmax.' endP=[(i-1),piancha,xdpiancha,xk,yk]';例2.6.7 判断x=0是方程f(x)=2e3x-9x2-6x-2=0的几重根？在区间*[0,1]上，分别用牛顿切线法和求重根的修正牛顿切线公式求此根的近似值xk，使精确到ε=10-4.解 经判断知，x=0是方程f(x)=0的三重根.在MATLAB工作窗口输入程序>> [k,piancha,xdpiancha,xk,yk]= newtonqx (0.5,1e-4, 1e-4,100)运行后整理结果列入表 2-20.* 迭代次数k=20，精确到ε=10的根的近似值是xk =0.000 2，其函数值是yk =5.755 8e-011，是一阶（线性）收敛速度.根据求重根的修正牛顿切线公式，在MATLAB工作窗口输入程序>> [k,piancha,xdpiancha,xk,yk]= newtonxz(3,0.5,1e-4, 1e-4,100)运行后整理结果得表 2-21. -47 0.000 10 迭代次数k=7，精确到ε=10的根的近似值是xk =1.121 1e-006，piancha = 9.975 3e-005，xdpiancha=88.974 8，其函数值是yk = -8.881 8e-016，是二阶收敛（平方收敛）速度.可见，求重根的修正牛顿切线公式比牛顿切线法收敛速度快得多.-488.974 84 0.000 00 -8.881 78e-016function [k,piancha,xdpiancha,xk,yk]=newtonxz1(x0,tol,ftol,gxmax) x(1)=x0;for i=1: gxmaxu(i)=fnq(x(i))/dfnq(x(i));du(i)=1-fnq(x(i))*ddfnq(x(i))/((dfnq(x(i)))^2+eps); x(i+1)=x(i)-u(i)/du(i); piancha=abs(x(i+1)-x(i));xdpiancha=piancha/( abs(x(i+1))+eps); i=i+1; xk=x(i);yk=fnq(x(i)); if ((pianchak=i-1; xk=x(i); [(i-1) piancha xdpiancha xk yk] end endif i>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax.') endP=[(i-1),piancha,xdpiancha,xk,yk]';例2.6.8 用未知重数的求重根的修正牛顿切线法，求方程 在区间[0,1]上的根，使精确到ε=10，并且与例2.6.7的结果比较. 解 在MATLAB工作窗口输入程序>> [k,piancha,xdpiancha,xk,yk]=newtonxz1(0.5,1e-4, 1e-4,100) 运行后整理结果得表 2-22.表 2-22-4piancha xdpiancha xk yk 0.599 23 6.038 57 -0.099 23 -0.008 18 0.096 98 43.054 79 -0.002 25 -0.000 00 0.002 25
-0.000 00 0.000 00 0.000 02 0.945 62 -0.000 03 -0.000 00 即迭代k =4次，piancha =2.461 55e-005， xdpiancha =0.945 62， xk =-2.603 09e-005， yk =-1.325 61e-013.这些结果与例2.6.7的结果比较见下表 2-23.k 1 2 3 4可见，未知重数的修正牛顿切线法的迭代次数和偏差最小；已知重数的修正牛顿切线法的迭代次数比牛顿切线法少13次，比未知重数的修正牛顿切线法多3次，但是根的近似值最接近根0；用牛顿切线法求得的近似根是三者中与根0距离最远的. 2.7 割线法及其2.7 割线法及其MATLAB程序 程序 割线法又称弦截法、弦位法、弦割法和弦法.牛顿切线法的突出优点是收敛速度快，但它也有明显的缺点：公式中含有导数，当f(x)较复杂时，使用不方便.为了避免切线法计算导数的麻烦，我们现在设法利用一些函数值f(xk)，f(xk-1)，…来回避导数值f'(xk)的计算.下面具体介绍割线法. 2.7.2 割线法的2.7.2 割线法的MATLAB程序 用割线法求方程f(x)=0的近似根x（精度为tol）需要自行编制程序. 01function [k,piancha,xdpiancha,xk,yk]=gexian (x01,x02,tol,ftol,gxmax) x(1)=x01;x(2)=x02; for i=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(i)/( v(i)); 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-2) piancha xdpiancha xk yk] if (abs(yk)if i>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax.') endP=[(i-2),piancha,xdpiancha,xk,yk]';2x2例2.7.1 用割线法求方程f(x)=e-3x=0的一个实根的近似值xk，使精确到.解 先用作图法确定初始值x01和x02，然后在MATLAB工作窗口输入程序>> [k,piancha,xdpiancha,xk,yk]= gexian (-0.4,-0.3,1e-4, 1e-4,100)运行后输出结果经整理得表 2-24.由此可见，迭代k =3次，得到精度为10-3的根的近似值xk =-0.390 6，其函数值为yk =3.860 7e-008，xk的相邻偏差为piancha =3.327 1e-005，其相误差为xdpiancha =8.516 8e-005.表 2-24k 1 2 3 piancha xdpiancha xk yk 0.090 090.230 95 -0.390 09 0.000.000 590.001 51 -0.390 68 -0.000 0.000 030.000 09 -0.390 65 0.000 00 2.8 抛物线法及其2.8 抛物线法及其MATLAB程序 程序 2.8.2 抛物线法的2.8.2 抛物线法的MATLAB程序 程序用抛物线法求方程f(x)=0的近似根x（精度为tol）需要自行编制程序.012的M文件function [k,piancha,xdpiancha,xk,yk]=paowu(x1,x2, x3,tol,ftol,gxmax)X=[x1, x2, x3]; for i=1: gxmaxh0= X(1)- X (3); h1= X (2)- X (3); u0= fnq (X (1))- fnq (X (3) ); u1= fnq (X (2))- fnq (X (3)); c= fnq (X (3));fenmu= h0^2* h1- h1^2* h0; afenzi= u0* h1- u1* h0; bfenzi= u1*h0^2-u0*h1^2;a= afenzi/ b= bfenzi/ deta=b^2-4*a*c;cha=[abs(X(3)-X(1)),abs(X(3)-X(2)),abs(X(2)- X(1))];piancha =min(cha); xdpiancha= piancha/( abs (X(3))+eps); if detadisp('提示：根的判别式值为负数.'),detakf=sqrt(deta); elseif detadisp('提示：根的判别式值为零.'),detakf=0; elsedisp('提示：根的判别式值为正数.'),detakf=sqrt(deta); end if bdetakf=-endz=-2*c/(b+ detakf);xk= X(3)+q1= xk - X (1); q2= xk - X (2); q3= xk - X (3); if abs(q2)X1=[X(2), X(1), X(3)]; X= X1;Y= fnq(X); endif abs(q3)X2=[X(1), X(3), X(2)]; X= X2;Y= fnq(X); endX(3)= Y(3)=fnq(X(3));yk= Y(3); [i piancha xdpiancha xk yk]if (abs(yk)if i>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值x0') endP=[i,piancha,xdpiancha,xk,yk]'; 例2.8.1 用抛物线法求方程f(x)=e2x-3x2=0的一个实根的近似值xk，使精确到ε=10-4.解 先用作图法确定初始值x01，x02和x03，然后在MATLAB工作窗口输入程序>> [k,piancha,xdpiancha,xk,yk]= paowu (-0.4,-0.3, -0.5,1e-4, 1e-4,100)运行后输出结果提示：根的判别式值为正数. ans =Columns 1 through 41.00 0.00 0.00 -0.39 Column 5-0.99提示：根的判别式值为正数. ans =Columns 1 through 42.00 0.61 0.38 -0.94 Column 5-0.32提示：根的判别式值为正数. ans =Columns 1 through 43.00 0.45 0.90 -0.59 Column 50.00 k = 3 piancha =1.676e-005 xdpiancha =4.760e-005 xk =-0.59 yk =2.313e-016即，迭代k =3次，得到精度为10-4的根的近似值xk =-0.390 6，其函数值为yk =2.220 4e-016，xk的相邻最小偏差为piancha =1.712e-005，其相对误差为xdpiancha =4.383 0e-005.将抛物线法与割线法的迭代结果进行比较，见表 2-25.显然，抛物线法比割线法的迭代 表 2-25迭代次抛物线法1. 66-0.000 06 2. 65-0.000 00 3. 650.000 00割线法-0.390 09 0.001 51 -0.390 650.001 81 -0.390 68 0.000 00 2.9 求解非线性方程组的牛顿法2.9 求解非线性方程组的牛顿法及其求解非线性方程组的牛顿法及其MATLAB程序 程序 求解非线性方程的牛顿法可以推广到二维非线性方程组的情形，也可以推广到更高维的情形. 2.9.1 求解2.9.1 求解二元求解二元非线性方程组的牛顿法二元非线性方程组的牛顿法及其非线性方程组的牛顿法及其MATLAB程序 程序用牛顿法求二元非线性方程组（2.40）的近似根xk,yk（精度为tol）需要自行编制程序.解二元非线性方程组的牛顿法二元非线性方程组的牛顿法的非线性方程组的牛顿法的MATLAB主程序 主程序输入的量：初始值X =[x0,y0],要求的近似根(xk,yk)的误差限tol， (xk,yk)的函数值fi(xk,yk),i=1,2的误差限ftol，迭代次数的最大值ngmax和函数fi(x,y),i=1,2及其一阶偏导数；输出的量：迭代ci=k次数得到的迭代向量值X=(xk,yk) 、向量的函数值Y={f1(xk,yk),f2(xk,yk)}、Y的2范数hanfan、迭代方程组（2.44）的系数行列式D的值、相邻两次迭代向量的2范数danfan=norm(Xk+1-Xk)和它的2 范数的相对误差xddf=danfan/norm(Xk+1)的值.(当迭代次数超过最大值ngmax时运行后输出信息’请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值x0’；当D=0时，运行后输出信息’请注意！迭代方程组（2.44）的系数行列式的值等于零.’)使用两个初始值x0,y0，根据公式（2.45），（2.46），（2.47）和（2.49），现提供名为newtonzu2.m的M文件function [ci,D,danfan,xddf,hanfan,Xk,Yk]=newtonzu2(X,tol,ftol,ngmax) Y=Z(X);for i=1:ngmaxdY=JZ(X);D=det(dY);A1=[Y(1),dY(1,2); Y(2),dY(2,2)]; A=det(A1);B1=[dY(1,1), Y(1); dY(2,1), Y(2)];B=det(B1);Xk=X-[A,B]/D; hanfan=norm(Y);danfan=norm(Xk-X); xddf=danfan/(norm(Xk)+eps); X=Xk;Y=Z(X);ci=i; if D~=0ci=i; Xk=X-[A,B]/D;Yk=Y; [ci,D,danfan, xddf, hanfan , X, Y]; elsedisp('请注意！迭代方程组（2.44）的系数行列式的值等于零.') endif (hanfan [ci,D,danfan, xddf, hanfan , X, Y]; endif i>gxmax'请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值x0.' end end22??x+y=16,例2.9.1 用迭代公式（2.46）求解方程组 ?2 2??x-y=2.解 这里介绍两种求解已知方程组的方法.方法1 用解1 用解二元用解二元非线性方程组的牛顿法二元非线性方程组的牛顿法的非线性方程组的牛顿法的MATLAB主程序求解主程序求解.求解. （1）设f1(x)=x+y-16，22f2(x)=x2-y2-2，则?f1(x,y)?x?y?x（2）作函数f1(x)和f2(x)的图形.输入程=2x，?f1(x,y)=2y，?f2(x,y)=2x，?f2(x,y)=-2y.?y序>>syms x yF1=x^2-y^2-2;F2=x^2+y^2-16; ezplot(F1,[-4.5,4.5]), hold on ezplot(F2,[-4.5,4.5]) ,hold off运行后屏幕显示图2-29.取两个初始值x0=2,y0=2.（3）建立并保存名为Z.m的M文件function F=Z(X)x=X(1);y=X(2); F=zeros(1,2);F(1)= x^2+y^2-16; F(2)= x^2-y^2-2; （4）建立并保存名为JZ.m的M文件 function dF=JZ(X)x=X(1);y=X(2); dF=zeros(2,2); dF(1,1)=2*x; dF(1,2)=2*y; dF(2,1)=2*x; dF(2,2)=-2*y;（5）保存名为newtonzu2.m的M文件； 图2-29 （6）在MATLAB工作窗口输入程序:>> X=[2,2]; [k,D,danfan, xddf, hanfan , Xk, Yk]= newtonzu2 (X,1e-4,1e-4,100) （7）运行后输出结果k = 5 D =-63.93 danfan =3.168e-011 xddf =9.920e-012 hanfan =3.849e-010 Xk =2.68 2.49 Yk =1.0e-015 *0 -0.13 方法2 用解矩阵方程的方法求解2 用解矩阵方程的方法求解.用解矩阵方程的方法求解. 编写如下程序while(norm(u)>tol*norm([x(i),y(i)]’))A=2*[x(i),y(i);x(i),-y(i)];b=[16-x(i)^2-y(i)^2,2-x(i)^2+y(i)^2]’;u=A\b;x(i+1)=x(i)+u(1); y(i+1)=y(i)+u(2); i=i+1;k(i)=i;if(i> gxmax)error('请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值');endend[k’,x’,y’]（2）运行后输出结果,取初始值x=(2,2)T，迭代次数n=6，精度ε=10-，6x6=3.=2.645751与精确解的误差已不超过ε. 2.9.2 求解2.9.2 求解n元非线性方程组的牛顿法及其非线性方程组的牛顿法及其MATLAB程序解元非线性方程组的牛顿法的非线性方程组的牛顿法的MATLAB主程序 主程序输入的量：初始值X0,要求的近似根x(k)的误差限tol， x(k)的函数值fi(x(k)),i=1,2,L,n的误差限ftol，迭代次数的最大值gxmax和函数fi(x(k)),i=1,2,L,n及其一阶偏导数；输出的量：迭代k次数得到的迭代向量值x(k)(k)、向量的函数值fi(x) i=1,2,L,n及其2范数hanfan、迭代方程组（2.51）的雅可比矩阵（2.52）的行列式D的值、相邻两次迭代向量的2 范数danfan=norm(x(k+1)-x(k))和它的2 范数的相对误差xddf= danfan/norm(x(k+1))的值.(当迭代次数超过最大值gxmax时运行后输出信息’请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值’；当D=0时，运行后输出信息’请注意！迭代方程组（2.51）的系数行列式的值等于零.’)使用两个初始值x0,y0，根据公式（2.50），（2.52）和（2.54），现提供名为newtonzu2.m的M文件function [ci,D,danfan, xddf, hanfan , Xk, Yk]= newtonzun (X, tol,ftol,gxmax)Y=Z(X);for i=1:gxmaxdY=JZ(X);D=det(dY);Xk=X-(dY\Y')';hanfan=norm(Y);danfan=norm(Xk-X);xddf=danfan/(norm(Xk)+eps);X=Xk;Y=Z(X);ci=i;if D~=0ci=i; Xk=X-(dY\Y')';Yk=Y;[ci,D,danfan,xddf,hanfan,X,Y];elsedisp('请注意！迭代方程组（2.51）的系数行列式的值等于零. ')endif (hanfan [ci,D,danfan, xddf, hanfan , X, Y];endendif i>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值.')end 29．?x+y=16,-6例2.9.2 用迭代公式（2.54）求解方程组 ?2要求精度. ε=102??x-y=2,解 这里介绍两种求解已知方程组的方法.方法1 用解1 用解n元非线性方程组的牛顿法的非线性方程组的牛顿法的MATLAB主程序求解主程序求解 求解在MATLAB工作窗口输入程序>> X=[2,2];[k,D,danfan, xddf, hanfan , Xk, Yk]= newtonzun (X,1e-6,1e-6,100)运行后输出结果k=5，D=-63.93，danfan=3.168e-011xddf=9.920e-012，hanfan=3.849e-010Xk = 3.00 2.59Yk = 1.0e-015 *0 -0.13方法2 用解矩阵方程的方法求解2 用解矩阵方程的方法求解在MATLAB工作窗口输入程序>> x0=2;y0=2; gxmax =10;tol=1e-6;x(1)=x0;y(1)=y0;i=1;u=[1 1];k(1)=1;while(norm(u)>tol*norm([x(i),y(i)]’))A=2*[x(i),y(i);x(i),-y(i)];b=[16-x(i)^2-y(i)^2,2-x(i)^2+y(i)^2]’;u=A\b; x(i+1)=x(i)+u(1); y(i+1)=y(i)+u(2); i=i+1;k(i)=i;if(i> gxmax)error('请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值’);endend[k’,x’,y’]运行后输出结果ans =1.00 2.00 2.002.00 3.00 2.003.00 3.38 2.274.00 3.32 2.805.00 3.32 2.696.00 3.00 2.596即，取初始值x=(2,2)T，迭代次数n=6，精度ε=10-，x6=3.=2.645751与精确解的误差已不超过ε. 2.9.3 求解2.9.3 求解n元非线性方程组的拟牛顿法及其非线性方程组的拟牛顿法及其MATLAB程序 用拟牛顿法求n元非线性方程组（2.50）的近似根x(k) （精度为tol），需要根据给定的方程组（2.50）和（2.58）分别自行编制名为NN.m和NF.m的M文件作为调用函数程序，然后根据（2.57）式，编写名为ninewton.m的主程序解元非线性方程组的拟牛顿法的非线性方程组的拟牛顿法的MATLAB程序 程序输入的量：初始值X01,X02，要求的近似根x(k)的误差限tol， x(k)的函数值fi(x(k)),i=1,2,L,n的误差限ftol，迭代次数的最大值gxmax和函数fi(x(k)),i=1,2,L,n及其一阶偏导数.),i=1,2,L,n及其2范数hanfan、迭代方程组（2.51）的雅可比矩阵（2.52）的行列式D的值、相邻两输出的量：迭代k次数得到的迭代向量值x 30． (k)、向量的函数值fi(x(k)次迭代向量的2范数danfan=norm(x(k+1)-x(k))和它的2范数的相对误差xddf=danfan/norm(x(k+1))的值.(当矩阵（2.58）行列式的值等于零时，输出信息’请注意！迭代方程组（2.59）中的矩阵（2.58）行列式的值等于零.’；当迭代次数超过最大值gxmax时运行后输出信息’请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值’) .function [k,hanfan,danfan,xddf,Xk,Yk]=ninewton(X01,X02,tol, ftol,gxmax)X=[X01;X02];n=length(X);for j=1:gxmaxY=NN(X); A=NF(X);D=det(A); Xk=X-((A+ones(n)*eps)\Y’)’;hanfan=norm(Y(2));danfan=norm(Xk(2)-X(2));xddf=danfan/(norm(Xk(2))+eps);j=j+1;X=Xk;Y=Z(X);Yk=Y;k=j;warning off MATLAB:singularMatrixif D~=0Xk=X-(A\Y’)’;Yk=Y;j=j+1; [j-1,D,danfan, xddf, hanfan ,Xk, Yk]’elsedisp(‘请注意！迭代方程组（2.59）中的矩阵（2.58）行列式的值等于零.’) endif (hanfan [j-1,D,danfan, xddf, hanfan ,Xk, Yk]’;warning off MATLAB:singularMatrixendendif(j> gxmax)error(‘请注意：迭代次数超过给定的最大值gxmax.’);end31．本文由(www.wenku1.com)首发，转载请保留网址和出处！
免费下载文档：}