MATLAB应用 求解非线性方程
第7章 求解非线性方程
7.1 多项式运算在MATLAB 中的实现
一、多项式的表达
n 次多项式表达为：n a +??++=x a x a x a p(x )1-n 1-n 1n 0，是n+1项之和 在MATLAB 中，n 次多项式可以用n 次多项式系数构成的长度为n+1的行向量表示
[a0, a1,……an-1,an]
二、多项式的加减运算 设
有
两
个
多
项
式
n
a +??++=x a x a x a p1(x )1-n 1-n 1n 0和
m b +??++=x b x b x b p2(x )1-m 1-m 1m 0。它们的加减运算实际上就是它们的对应系
数的加减运算。当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。
例2 计算()()1635223-+++-x x x x a=[1, -2, 5, 3]; b=[0, 0, 6, -1]; c=a+b
例3 设()6572532345++-+-=x x x x x x f ，()3532-+=x x x g ，求f(x)+g(x)
f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; g1=[0, 0, 0, g];%为了和f 的次数找齐
f+g1, f-g1
三、多项式的乘法运算 conv(p1,p2)
例4 在上例中，求f(x)*g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; conv(f, g)
四、多项式的除法运算
[Q, r]=deconv(p1, p2)
表示p1除以p2，给出商式Q(x)，余式r(x)。Q,和r 仍为多项式系数向量 例4 在上例中，求f(x)/g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; [Q, r]=deconv(f, g) 五、多项式的导函数
p=polyder(P)：求多项式P 的导函数 p=polyder(P,Q)：求P·Q 的导函数
[p,q]=polyder(P,Q)：求P/Q 的导函数，导函数的分子存入p ，分母存入q 。 参数P,Q 是多项式的向量表示，p,q 也是多项式的向量表示。
例4 求有理分式()100
765105
853*********--+++-+-+=x x x x x x x x x x f 的导函数
P=[3, 5, 0, -8, 1, -5]; %有理分式分子 Q=[10, 5, 0, 0, 6, 0, 0, 7, -1, 0, -100]; %有理分式分母 [p,q]=polyder(P,Q) 六、多项式求根
多项式求根就是求满足多项式p(x)＝0的x 值。N 次多项式应该有n 个根。这些根可能是实根，也可能是若干对共轭复根。其调用格式是
x=roots(P)
其中P 为多项式的系数向量，求得的根赋给向量x ，即x(1),x(2),…,x(n)分别代表多项式的n 个根。
该命令每次只能求一个一元多项式的根，该指令不能用于求方程组的解，必须把多项式方程变成P n (x) = 0的形式；
例4 求方程123+=x x 的解。
首先将方程变成P n (x) = 0的形式：0123=--x x
roots([1 -1 0 -1])
例5 求多项式x 4+8x 3-10的根。 A=[1,8,0,0,-10]; x=roots(A)
若已知多项式的全部根，则可以用poly 函数建立起该多项式，其调用格式为：
P=poly(x)
若x 为具有n 个元素的向量，则poly(x)建立以x 为其根的多项式，且将该多项式的系数赋给向量P 。
例6 已知 f(x)=3x 5+4x 3-5x 2-7.2x+5 (1) 计算f(x)=0 的全部根。
(2) 由方程f(x)=0的根构造一个多项式g(x)，并与f(x)进行对比。 P=[3,0,4,-5,-7.2,5];
X=roots(P) %求方程f(x)=0的根 G=poly(X) %求多项式g(x) 将这个结果乘以3，就与f(x)一致
7.2 求解非线性方程f ( x ) = 0
方程求根的一般形式是求下列方程的根：
f ( x ) = 0 (l)
实际上，就是寻找使函数 f ( x ）等于零的变量x ，所以求方程（l ）的根，
也叫求函数 f ( x ）的零点。如果变量x 是列阵，则方程（l ）就代表方程组。
当方程（l ）中的函数 f (x ）是有限个指数、对数、三角、反三角或幂函数的组合时，则方程（l ）被称为超越方程，例如 e -x - sin （πx / 2 ) ＋lnx = 0 就是超越方程。
当方程（l ）中的函数f （x ）是多项式时，即 f （x ）＝P n （x ）= a n x n + a n-1x n ＋ … ＋ a l x + a 0，则方程（l ）就成为下面的多项式方程，也称代数方程：
P n （x ）= a n x n + a n-1x n ＋ … ＋ a l x + a 0 = 0 ( 2 )
P n （x ）的最高次数n 等于2、3时，用代数方法可以求出方程（2）的解析解，但是，当n ≥ 5时，伽罗瓦（Galois ）定理已经证明它是没有代数求根方法的。至于超越方程，通常很难求出其解析解。所以，方程（l ）的求解经常使用作图法或数值法，而计算机的发展和普及又为这些方法提供了广阔的发展前景，使之成为科学和工程中最实用的方法之一。
本章首先介绍求解 f ( x ) = 0 的 MATLAB 符号法指令，然后介绍求方程数值解的基本原理，最后再介绍求解 f ( x ) = 0 的 MATLAB 数值法指令。
一、符号方程求解
在MATLAB 中，求解用符号表达式表示的代数方程可由函数solve 实现，其调用格式为：
solve(s)：求解符号表达式s 的代数方程，求解变量为默认变量。 当方程右端为0时，方程可以不标出等号和0，仅标出方程的左端。 solve(s,v)：求解符号表达式s 的代数方程，求解变量为v 。
solve(s1,s2,…,sn,v1,v2,…,vn)：求解符号表达式s1,s2,…,sn 组成的代数方程组，求解变量分别v1,v2,…,vn 。
例1． 解下列方程。
1．2
2144212-+
=-++x x x x x= solve('1/(x+2)+4*x/(x^2-4)=1+2/(x-2)', 'x')
2．17433=---x x x f=sym('x-(x^3-4*x-7)^(1/3)=1') x= solve(f)
3．143sin 2=??? ??
-πx
x= solve('2*sin(3*x-pi/4)=1') 4．010=-+x xe x
x= solve('x+x*exp(x)-10', 'x') %仅标出方程的左端 二、求方程f ( x ) = 0数值解的基本方法
并非所有的方程 f ( x ) = 0 都能求出精确解或解析解，不存在这种解的方程
就需要用数值解法求出近似解，有几种常见的数值解法基本原理：二分法。
1 求实根的二分法原理
设方程 f (x) =0中的函数 f ( x ）为实函数，且满足：
① 函数 f (x ）在[ a , b]上单调、连续；
② 方程 f (x) = 0 在（a , b ）内只有一个实根 x*。
则求方程 f (x) = 0 的根，就是在（a, b ）内找出使f (x ）为零的点x*：f (x*) = 0 ，即求函数 f ( x ) 的零点。因为 f (x ）单调连续，由连续函数的性质可知，若任意两点aj ，bj ∈[ a , b] ，而且满足条件 f (aj) f (bj) < 0 ，则闭区间[aj , bj] 上必然存在方程的根x*，即 x*∈[aj , bj]。
据此原理提出求实根的二分法如下图所示，
图1 方程求根二分法原理示意图
先用中点2
1b
a b +=
将区间[a, b]平分为两个子区间 (a,b 1）和（b 1, b)，方程的根必然在子区间两端点上函数值之积小于零的那一半中，即不在(a ,b 1）内，就在（b 1 ,b )内，除非 f(b 1) = 0 ，于是寻根的范围缩小了一半。图1中的根x*在区间中点左侧，即 x*∈(a , b l ）。再将新的含根区间( a , b 1）分成两半，重复上述步骤确定出更新的含根子区间。如此重复n 次，设含根区间缩小为（a n , b n ），则方程的根x*∈（a n , b n ）, 这一系列含根的子区间满足：
( a , b ) D ? ( a l , b l ) ? ( a 2 , b 2 ) ? … ? ( a 0， b 0）? …
由于含根区间范围每次减半，子区间的宽度为n n n a
b a b 2
-=- (n = 1,2,….)，显然
当n →∞时，(b n 一a n ）→0，即子区间收敛于一点x*，这个点就是方程的根。若n 为有限整数，取最后一个子区间的中点2
n
n n b a x +=
作为方程根的近似值，它
满足 f ( x n ）≈0 ，于是有：
12
221*+-=-?≤
-n n n a b a b x x 这就是近似值x n 的绝对误差限。假定预先要求的误差为ε，由1
2
+-<
n a
b ε便可以求出满足误差要求的最小等分次数n 。 下面是二分法的程序
function [c,err,yc] =bisect (f,a,b,delta)
%Input - f is the function input as a string ‘f ’ % - a and b are the left and right end points %. - delta is the tolerance %Output - c is the zero % - yc=f(c)
% - err is the error estimate for c ya=feval (f,a); yb=feval (f,b);
if ya*yb>0, break, end %表示无解，结束
maxl=l+round( (log (b-a) -log (delta))/log (2)); %从误差表达式得到最小等分次数n for k=1:max1
c=(a+b)/2; %取区间中点 yc=feval (f,c); if yc==0 a=c;
b=c; %这时解已经找到 elseif yb*yc>0 b=c; %区间减半 yb=yc; else ?a=c; ya=yc; end
if b-a < delta, break, end end
c=(a+b)/2; err=abs(b-a); yc=feval (f, c)
2 迭代法
迭代法是计算数学中的一种重要方法，用途很广，求解线性方程组和矩阵特征值时也要用到它。这里结合非线性方程的迭代法求解，介绍一下它的基本原理。
迭代法基本原理
迭代法的基本原理就是构造一个迭代公式，反复用它得出一个逐次逼近方程根的数列，数列中每个元素都是方程根的近似值，只是精度不同。
迭代法求解方程
f ( x ) = 0 （1） 时，先把方程等价地变换成形式
f ( x ) = x －g(x) = 0 ， （2） 移项得出：
x = g(x) ( 3 )
若函数g (x ）连续，则称（3）为迭代函数。用它构造出迭代公式： x k+1= g ( x k ） , k = 0 , l , 2 ， … ( 4 ) 从初始值 x 0出发，便可得出迭代序列：
{ x k ｝= x 0, x 1, x 2,….x k ,….. （ 5 ） 如果迭代序列（5 ）收敛，且收敛于x*，则由式（4）有：
()()()()()0***lim 1==-=-+∞
→x f x x g x x g k k
可见 x*便是方程（l ）的根。 迭代法几何意义：
如下图所示，解方程f ( x ) = 0可以等价地变换成求解 x = g ( x )，
图 4 - 2 方程求根迭代法原理示意图
在几何上，就等价求曲线y ＝x 和y ＝g ( x ）交点P*的坐标 x*。求迭代序列（5) ，就等于从图中x 0点出发，由函数y ＝g ( x 0）得出y ＝P 0，代入函数y ＝x 中得出Q 1，再把Q 1的x 坐标 x 1代入方程y= g ( x ）得出P 1，如此继续下去，便可在曲线y ＝g ( x )上得到一系列的点P 0，P 1， … ，P k ， … ，这些点的x 坐标便是迭代数列 x l , x2 ， … ， x k ， … ，它趋向于方程（l ) 的根 x*，数列的元素就是方程根的近似值。数列的收敛就等价于曲线y ＝x 和y ＝g ( x ）能够相交于一点。
迭代公式收敛定理
要想用迭代法求出方程根的近似值，迭代序列（4 - 5）必须收敛。下面的定理给出了迭代法的收敛条件，同时也给出了迭代公式的误差。
收敛定理：方程 x = g ( x ）在（ a , b ）内有根 x*，如果： ① 当x ∈[a,b]时，g( x ）∈[a,b];
② g ( x ）可导，且存在正数 q < 1，使得对于任意x ∈[a,b]都有|g ’( x )|≤ q < 1，则有以下结论。
① 方程x = g ( x ）在（a , b ）内有唯一的根x*。
② 迭代公式 x k+1 = g ( x k ）对（a , b ）内任意初始近似根 x 0均收敛于x*。 ③ 近似根 x k 的误差估计公式为：
011*x x q
q x x k
n --≤- (4 - 6)
3 切线法
切线法就是从函数曲线上的一点出发，不断用曲线的切线代替曲线，求得收敛于根的数列。
切线法原理：
解非线性方程f( x ) = 0 的切线法也称牛顿法，它是把方程线性化的一种近似方法，用函数 f (x ）的切线代替曲线产生一个收敛于方程根的迭代序列，从而得到方程的近似根。
把函数 f ( x ）在某一初始值 x 。点附近展开成泰勒级数：
()()()()()
()
+''-+'-+=!
202
0000x f x x x f x x x f x f (4 - 7) 取其线性部分，近似地代替函数f (x ）可得方程的近似式：
()()()()0000='-+≈x f x x x f x f
设()00≠'x f ，解该近似方程可得：
()
()
0001x f x f x x '-
= 把函数 f ( x ）在x l 点附近展开成泰勒级数，取其线性部分替代函数f (x)，设()01≠'x f ，得：
()
()
1112x f x f x x '-
= 如此继续做下去，就可以得到牛顿迭代公式：
()
()
k k k k x f x f x x '-
=+1 ，k = 0 , l , 2 ， … ( 8 ) 由式（8）得出的迭代序列 x l , x 2 ， … ，x k … ，在一定的条件下收敛于方程的根x*。 2 ．几何意义
图 4 一 3 方程求根切线法原理示意图
选取初值x 0后，过()()00,x f x 点作曲线()x f y =的切线，其方程为
()()()000x f x x x f y '-=-。设切线与X 釉的交点为x 1，则()
()
0001x f x f x x '-
=，再过()()11,x f x 作切线，与x 轴的交点为()
()
1112
x f x f x x '-
=， 如此不断作切线，求与x 轴的交点，便可得出的一系列的交点x 1，x 2，…，x k ，…，它们逐渐逼近方程的根x*。
3 ．切线法的收敛性
理论可以证明，在有根区间[a, b]上，如果()00≠'x f 、()00≠''x f 连续且不变号，则只要选取的初始近似根x 0满足()()000>'''x f x f f ( x 。），切线法必定收敛。它的收敛速度经推导可得出：
()()
()21**2**x x x f x f x x k k -'''=
-+ ( 9)
()()
*2*x f x f '''是个常数，式(9)表明用牛顿迭代公式在某次算得的误差，与上次误差的平方成正比，可见牛顿迭代公式的收敛速度很快。 4 . 2 . 4 割线法（弦截法）
应用切线法的牛顿迭代公式时，每次都得计算导数()k x f '，若将该导数用差商代替，就成为割线法（有时称快速弦截法）的迭代公式：
()()
()k k k k k k k x f x f x f x x x x 11
1--+---
=，k = 0 , l , 2 ， … ( 4 一 10 )
割线法的几何意义也很明显。如图 所示，
p 0
p 1
p 2
p 3x 0
x 1x 2x 3
x 4p 4
图 4 一 4 方程求根割线法原理示意图
过点（x 0，f( x 0））和（x 1，f (x 1)）作函数y ＝f(x ）曲线的割线，交X 轴于点x 2，再过点（x 1，f (x 1)）和（x 2，f (x 2)）作曲线的割线，交X 轴于点x 3，……一直做下去，则得到一系列割线与 X 轴的交点，这些交点序列将趋于方程的根 x*。
非线性方程的数值解法还有许多，这里仅介绍了几种基本方法的原理。 二分法简单方便，但收敛速度慢；
迭代法虽然收敛速度稍微快点，但需要判断能否收敛； 只要初值选取得当，切线法具有恒收敛且收敛速度快的优点，但需要求出函数的导数；
弦截法不需要求导数，特别是前面介绍的快速弦截法，收敛速度很快，但是需要知道两个近似的初始根值才能作出弦，要求的初始条件较多。
这些方法各有千秋，需根据具体情况选用。
三、方程f(x) = 0数值解的MATLAB实现
MATLAB中求方程数值解的办法很多，有的是专用指令，有的是根据方程性质而借用其他专用指令求得的。
4 . 3 . 2 求函数零点指令fzero
求解方程f ( x ) = 0的实数根也就是求函数f ( x）的零点。MATLAB中设有求函数f (x）零点的指令fzero，可用它来求方程的实数根。该指令的使用格式为：fzero (fun, x0, options)
①输入参数fun为函数f (x）的字符表达式、内联函数名或M函数文件名。
②输入参数x0为函数某个零点的大概位置（不要取零）或存在的区间［x i，x j］，要求函数f (x）在x0点左右变号，即f (x i)f (x j) < 0。
③输入参数options可有多种选择，若用optimset ('disp', 'iter'）代替options 时，将输出寻找零点的中间数据。
④该指令无论对多项式函数还是超越函数都可以使用，但是每次只能求出函数的一个零点，因此在使用前需摸清函数零点数目和存在的大体范围。为此，一般先用绘图指令plot, fplot或ezplot 画出函数f (x）的曲线，从图上估计出函数零点的位置。
例 4 一 5 求方程x2 + 4sin(x) = 25 的实数根（-2π＜x < 2π）。
解：
一、fun为函数f (x）的字符表达式
（l）首先要确定方程实数根存在的大致范围。为此，先将方程变成标准形式f(x) = x2 + 4sin(x) - 25 = 0。作f(x)的曲线图：
x=-2*pi:0.1:2*pi;
f=x.^2+4*sin(x)-25;
plot(x,f);grid on;
从曲线上可以看出，函数的零点大约在x1≈ - 4和x2≈ 5附近。
(2）直接使用指令fzero 求出方程在x1≈ - 4时的根。
x1= fzero ('x^2+4*sin(x)-25',-4)
若键入：fzero ('x^2+4*sin(x)-25',-4, optimset('disp', 'iter'))，将显示迭代过程。
中间数据表明，求根过程中不断缩小探测范围，最后得出- 4附近满足精度的近似根。
(3）求x2≈ 5的根：
x2= fzero ('x^2+4*sin(x)-25',5)
二、fun为函数f (x）的M函数文件名
将方程x2 + 4sin(x) = 25编成M函数文件（实用中在函数较为复杂、而又多次重复调用时，才这样做），用fzero求解。
(1)在M文件编辑调试窗中键入：
function yy = li4_5(x)
yy= x^2+4*sin(x)-25;
以li4_5为文件名存盘，退出编辑调试窗，回到指令窗。
(2）确定根的大体位置;
(3) 在指令窗中键入下述指令可求出- 4 附近的根：
x1= fzero (' li4_5',-4)
键入下述指令可求出5附近的根：
x2= fzero (' li4_5',5)
三、fun 为函数f (x ）的内联函数名
内联函数是MATLAB 提供的一个对象(Object)。它的性状表现和函数文件一样，但内联函数的创建比较容易。
inline('CE')
'CE'是字符串，CE 为不包含赋值符号“=”的表达式。 上式把串表达式转化为输入宗量自动生成的内联函数。 上述调用格式将自动地对CE 进行辫识，把CE 中由字母／数字组成的连续字符认做变量，除“预定义变量名（如 i , j , pi ) ”和“常用函数名（如 sin )”以外的由字母／数字组成的连续字符将被认做变量。但注意：若连续字符后紧接“左圆括号”，那么将不被当作输入宗量。如 x ( 1 ) ，就不会认做输入宗量处理。
inline('CE',arg1,arg2,…)
上述调用格式把串表达式转化为arg1,arg2等指定输入宗量的内联函数；这种调用格式是创建内联函数的最稳妥、可靠途径。输入宗量字符可表达得更自如。
将函数f (x ）写成内联函数的形式： f = inline ('x^2+4*sin(x)-25') 这时内联函数名为f
分别求x 1 ≈ - 4和x 2 ≈ 5时的根： x1= fzero (f,-4), x2= fzero (f,5)
例4-6 求f(x)=x-10x +2=0在x0=0.5附近的根。
从f(x)的曲线看(x=-2.5:0.01:0.5; fx=x-10.^x+2; plot(x,fx))，
曲线的零点有两个，一个在x=-2附近，另一个在x=0.5附近
(1) 建立函数文件funx.m 。(function [输出变量列表]=函数名(输入变量列表)) function fx=funx(x) fx=x-10.^x+2;
(2) 调用fzero 函数求根。 z=fzero('funx',0.5)
例2 求()51
+-=x
x x f 在50-=x 和10=x 作为迭代初值时的零点。
从f(x)的曲线看(x=-7:0.01:2; f= x-1./x+5; plot(x,f))，
曲线的零点有两个，一个在x=-5附近，另一个在x=1附近 (1) 建立函数文件fz.m 。 function f=fz(x) f= x-1./x+5;
(2) 调用fz.函数求根。
fzero('fz',-5) %以-5作为迭代初值 fzero('fz',1) %以1作为迭代初值
7. 3 求解非线性方程组数值解的迭代法
一、符号方程组求解
在MATLAB 中，求解用符号表达式表示的方程组仍然可由函数solve 实现，其调用格式与解用符号表达式表示的方程一样。
例 解下列方程组。
1．???
????=+=+41128113
3y x y x
[x y]=solve('1/x^3+1/y^3=28', '1/x+1/y=4', 'x,y')
2．???=+=+2
9833v u v u
[u v]=solve('u^3+v^3=98', 'u+v=2', 'u,v')
3．?
??=+=+29833y x y x [x y]=solve('x+y=98', 'x^(1/3)+y^(1/3)=2', 'x,y') 回车后出现下面的提示
Warning: Explicit solution could not be found.
> In D:\MATLAB6p5\toolbox\symbolic\solve.m at line 136
如果做代换：33,v y u x →→，方程3就变成方程2，就可解
这个问题说明，符号求解并不是万能的。
如果用MATLAB 得出无解或未找到所期望的解时，应该用其它方法试探求解。
4．?
??=--=+023252
222y xy x y x [x y]=solve('x^2+y^2=5', '2*x^2-3*x*y-2*y^2') %变量由默认规则确定 二、求解非线性方程组的基本方法
对于非线性方程组（以二元方程组为例，其他可以类推）
()()??
?==0,0
,2
1y x f y x f ( 11) 的数值解求法，跟一元非线性方程的切线法（牛顿法）雷同，也是把非线性函数
线性化，近似替代原方程得出数值解，所以也叫作牛顿迭代法。
假设方程组（11）的初始估计值为（x 0，y 0），可以把方程组（11）中的两个函数f 1(x ，y)和f 2(x ，y)在（x 0，y 0）处用二元泰勒级数展开，只取线性部分，移项得出：
()()()()()()()()()()???????-=-??+-??-=-??+-??0020002000200100010001,,,,,,y x f y y y y x f x x x y x f y x f y y y
y x f x x x y x f (12) 若系数矩阵行列式()
()()()0,,,,0020020010010≠????????=
y
y x f x
y x f y
y x f x
y x f J ，则方程组（12）的解为：
()()()
()
002002001001001,,,,1y x f y
y x f y x f x y x f J x x ????+
=，()
()
()()y
y x f y x f x
y x f y x f J y y ????+
=0020
020********,,,,1
方程组（11）中的两个函数f 1(x ，y)和f 2(x ，y)在(x 1，y 1)处，再用二元泰勒级数展开，只取线性部分，… …如此继续替代下去，直到方程组的根达到所要求的精度，就完成了方程组的求解。
求解方程组（11）还有许多其他办法，如“最速下降法”，它是利用方程组（11）构成所谓模函数()()[]()[]2
22
1,,y x f y x f x +=Φ，通过求模函数极小值的方法
得到方程组的数值解，诸如此类在此不再一一列举。
三、求方程组数值解的指令
fsolve 是用最小二乘法求解非线性方程组 F (X ）＝ 0 的指令，变量X 可
以是向量或矩阵，方程组可以由代数方程或者超越方程构成。它的使用格式为：
fsolve ( 'fun' , X0 , OPTIONS )
① 参数 fun 是编辑并存盘的M 函数文件的名称，可以用＠代替单引号对它进行标识。 M 函数文件主要内容是方程 F ( X ) = 0 中的函数 F (X) ，即方程左边的函数。
② 参数X0是向量或矩阵，为探索方程组解的起始点。求解将从X0出发，逐渐趋向，最终得到满足精度要求、最接近X0的近似根X* : F ( X*）≈ 0 。由于X0是向量或矩阵，无法用画图方法进行估计，实际问题中常常是根据专业知识、物理意义等进行估计。
③ 该指令输出一个与X0同维的向量或矩阵，为方程组的近似数值解。 ④ 参数OPTIONS 为设置选项，用它可以设置过程显示与否、误差、算法… …，具体内容可用 help 查阅。通常可以省略该项内容。
例 求方程组??
?
??=++==+++y z y x z xy x z y x 32107222222在x 0＝1 ，y 0＝1，z 0＝1附近的数值解。
解
（l ）在文本编辑调试窗中编辑 M 函数文件。首先将方程组变换成 F ( X )
= 0 的形式，??
?
??=+-+=-=+++-03020710222222z y y x z xy z y x x x, y, z 看成向量X 的三个分量。
function ms=li4_7(X)
ms(1) = X(1).^2 - 10* X(1) + X(2).^2 + X(3) +7; ms(2) = X(1).*X(2).^2 - 2* X(3);
ms(3) = X(1).^2 + X(2).^2 - 3* X(2) + X(3).^2;
这里，X (l ）＝x , X (2）＝y , X (3）＝z ，输出参量 ms 也有三个分量 用“li4_7”为M 函数文件名存盘。 (2）在指令窗中键入：
fsolve (' li4_7',[1 1 1]) 若键入：
x = fsolve (@li4_7,[1 1 1], optimset('Display', 'iter')) 则得出求解过程。
该方程也可以用 MATLAB 的符号指令solve 求解，但结果非常冗长。
例 求解方程组()()???
????
=++=+++-+=-131020006.1sin 1.08125.0cos 32
2πz e z y x yz x xy ，在 x 0＝0.1，y 0＝0.1和z 0
＝ - 0.1附近的数值解。
解：
首先将方程组变换成f j (x, y, z ）＝f (X ）＝0 ( j ＝1 , 2 , 3）的形式，设X 为一个三维向量，令X (l ）＝x , X (2）＝y , X (3）＝z ，则三维向量yy3＝f (X ）＝f j (x, y, z ），然后编程计算。
(l) 在文本编辑调试窗中编辑M 函数文件 function yy3=li4_8(X)
yy3(1) = 3*X(1) - cos(X(2)*X(3)) - 0.5;
yy3(2) = 2*X(1)^2 - 81*(X(2) + 0.1)^2 +sin(X(3)) + 0.06; yy3(3) = exp(-X(1)*X(2)) + 20*X(3) + 10*pi/3 -1; (2）在指令窗中键入： fsolve('li4_8',[0.1 0.1 -0.1])
这个方程组用符号指令solve 无法得出最终结果。
例1 求解方程组??
?
??==++=++000sin 2xyz z y x e z y x x 在(1,1,1)附近的解，并对结果进行验证。
(1) 建立函数文件myfun.m 。
function F=myfun(X) %X 、F 是性质相同的向量 x=X(1); y=X(2); z=X(3);
F(1)= sin(x)+y+z^2.*exp(x); F(2)=x+y+z; F(3)=x.*y.*z;
(2) 在给定的初值x0=1,y0=1,z0=1下，调用fsolve 函数求方程的根。 X=fsolve('myfun',[1,1,1],optimset('Display','off'))
(3) 将求得的解代回原方程，可以检验结果是否正确，命令如下： q=myfun(X)
例2 求非线性方程组???=+-=--0sin 3.0cos 6.00
cos 3.0sin 6.0y x y x x x 在(0.5,0.5) 附近的数值解。
(1) 建立函数文件myfun1.m 。
function F=myfun1(X) x=X(1); y=X(2);
F(1)=x-0.6*sin(x)-0.3*cos(y); F(2)=y-0.6*cos(x)+0.3*sin(y); function F=myfun1(X) x=X(1); y=X(2);
F(1)=x-0.6*sin(x)-0.3*cos(y); F(2)=y-0.6*cos(x)+0.3*sin(y);
(2) 在给定的初值x0=0.5,y0=0.5下，调用fsolve 函数求方程的根。 X=fsolve('myfun1',[0.5,0.5]',optimset('Display','off'))
将求得的解代回原方程，可以检验结果是否正确，命令如下： q=myfun1(X)
推荐-Broyden方法求解非线性方程组的Matlab实现 精品
Broyden方法求解非线性方程组的Matlab实现 注：matlab代码来自网络，仅供学习参考。 1.把以下代码复制在一个.m文件上 function [sol, it_hist, ierr] = brsola(x,f,tol, parms) % Broyden's Method solver, globally convergent % solver for f(x) = 0, Armijo rule, one vector storage % % This code es with no guarantee or warranty of any kind. % % function [sol, it_hist, ierr] = brsola(x,f,tol,parms) % % inputs: % initial iterate = x % function = f % tol = [atol, rtol] relative/absolute % error tolerances for the nonlinear iteration % parms = [maxit, maxdim] % maxit = maxmium number of nonlinear iterations % default = 40 % maxdim = maximum number of Broyden iterations % before restart, so maxdim-1 vectors are % stored % default = 40 % % output: % sol = solution % it_hist(maxit,3) = scaled l2 norms of nonlinear residuals % for the iteration, number function evaluations, % and number of steplength reductions % ierr = 0 upon successful termination % ierr = 1 if after maxit iterations % the termination criterion is not satsified. % ierr = 2 failure in the line search. The iteration % is terminated if too many steplength reductions % are taken. % % % internal parameter: % debug = turns on/off iteration statistics display as % the iteration progresses
MATLAB代码 解线性方程组的迭代法
解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;
if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵
MatLab求解线性方程组
MatLab解线性方程组一文通 当齐次线性方程AX=0,rank(A)=r 在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法. 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB 求解线性方程组 solve，linsolve 例： A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1]; %矩阵的行之间用分号隔开，元素之间用逗号或空格 B=[3;1;1;0] X=zeros(4,1);%建立一个4元列向量 X=linsolve(A,B) diff（fun，var，n）：对表达式fun中的变量var求n阶导数。 例如：F=sym（'u(x,y)*v(x,y)'）; %sym（）用来定义一个符号表达式 diff(F); %matlab区分大小写 pretty(ans) %pretty（）：用习惯书写方式显示变量；ans是答案表达式 非线性方程求解 fsolve(fun,x0,options) 为待解方程或方程组的文件名；fun其中 x0位求解方程的初始向量或矩阵； option为设置命令参数 建立文件fun.m： function y=fun(x) y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ... x(2) - 0.5*cos(x(1))+0.3*sin(x(2))]; >>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve')) 注： ...为续行符 m文件必须以function为文件头，调用符为@；文件名必须与定义的函数名相同；fsolve（）主要求解复杂非线性方程和方程组，求解过程是一个逼近过程。Matlab求解线性方程组 AX=B或XA=B 在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“/”和“\”。如： X=A\B表示求矩阵方程AX＝B的解； 的解。XA=B表示矩阵方程B/A＝X． 对方程组X＝A\B，要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数，方程X＝B/A同理。 如果矩阵A不是方阵，其维数是m×n，则有： m＝n 恰定方程，求解精确解； m>n 超定方程，寻求最小二乘解； m 中南大学MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券的自行百度。所需m文件照本文档做即可，即新建（FILE）→脚本（NEW-Sscript）→复制本文档代码→运行（会跳出保存界面，文件名默认不要修改，保存）→结果。第一题需要把数据文本文档和m文件放在一起。全部测试无误，放心使用。本文档针对做牛顿法求非线性函数题目的同学，当然第一题都一样，所有人都可以用。←记得删掉这段话 班级： ？ 学号： 姓名： 一、《MATLAB程序设计实践》Matlab基础 表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散 空间函数值来表示取向分布函数，是三维取向分布函数的一个实例。 由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一 个matlab程序画出如下的几种图形来分析其取向分布特征： （1）用Slice函数给出其整体分布特征； " ~ （2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25, 30, 35 … 90)切面上f分布情况(需要用到subplot函数)； (3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。 ( 备注：数据格式说明 解： （1）（ （2）将文件内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如 下： fid=fopen(''); for i=1:18 tline=fgetl(fid); end phi1=1;phi=1;phi2=1;line=0; f=zeros(19,19,19); [ while ~feof(fid) tline=fgetl(fid); data=str2num(tline); line=line+1;数据说明部分，与 作图无关此方向表示f随着 φ1从0,5,10,15, 20 …到90的变化而 变化 此方向表示f随着φ 从0,5,10,15, 20 … 到90的变化而变化 表示以下数据为φ2＝0的数据，即f（φ1，φ，0） 《MATLAB语言》课成论文 利用MATLAB求线性方程组 姓名：郭亚兰 学号：12010245331 专业：通信工程 班级：2010级通信工程一班 指导老师：汤全武 学院：物电学院 完成日期：2011年12月17日 利用MATLAB求解线性方程组 （郭亚兰 12010245331 2010 级通信一班） 【摘要】在高等数学及线性代数中涉及许多的数值问题，未知数的求解，微积分，不定积分，线性方程组的求解等对其手工求解都是比较复杂，而MATLAB语言正是处理线性方程组的求解的很好工具。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是讨论矩阵理论、与矩阵结合的有限维向量空间及其线性变换理论的一门学科。 【关键字】线性代数MATLAB语言秩矩阵解 一、基本概念 1、N级行列式A：A等于所有取自不同性不同列的n个元素的积的代数和。 2、矩阵B：矩阵的概念是很直观的，可以说是一张表。 3、线性无关：一向量组（a1,a2,…,an）不线性相关，既没有不全为零的数 k1,k2,………kn使得：k1*a1+k2*a2+………+kn*an=0 4、秩：向量组的极在线性无关组所含向量的个数成为这个向量组的秩。 5、矩阵B的秩：行秩，指矩阵的行向量组的秩；列秩类似。记：R(B) 遗传算法解非线性方程组的Matlab程序 程序用MATLAB语言编写。之所以选择MATLB，是因为它简单，但又功能强大。写1行MATLAB程序，相当于写10行C++程序。在编写算法阶段，最好用MATLAB语言，算法验证以后，要进入工程阶段，再把它翻译成C++语言。 本程序的算法很简单，只具有示意性，不能用于实战。 非线性方程组的实例在函数(2)nonLinearSumError1(x)中，你可以用这个实例做样子构造你自己待解的非线性方程组。 %注意：标准遗传算法的一个重要概念是，染色体是可能解的2进制顺序号，由这个序号在可能解的集合(解空间)中找到可能解 %程序的流程如下： %程序初始化，随机生成一组可能解(第一批染色体) %1: 由可能解的序号寻找解本身(关键步骤) %2：把解代入非线性方程计算误差，如果误差符合要求，停止计算 %3：选择最好解对应的最优染色体 %4：保留每次迭代产生的最好的染色体，以防最好染色体丢失 %5: 把保留的最好的染色体holdBestChromosome加入到染色体群中 %6: 为每一条染色体(即可能解的序号)定义一个概率(关键步骤) %7：按照概率筛选染色体(关键步骤) %8：染色体杂交(关键步骤) %9：变异 %10:到1 %这是遗传算法的主程序，它需要调用的函数如下。 %由染色体(可能解的2进制)顺序号找到可能解： %(1)x=chromosome_x(fatherChromosomeGroup,oneDimensionSet,solutionSum); %把解代入非线性方程组计算误差函数：(2)functionError=nonLinearSumError1(x); %判定程是否得解函数：(3)[solution,isTrue]=isSolution(x,funtionError,solutionSumError); %选择最优染色体函数： %(4)[bestChromosome,leastFunctionError]=best_worstChromosome(fatherChromosomeGroup,functionError); %误差比较函数：从两个染色体中，选出误差较小的染色体 %(5)[holdBestChromosome,holdLeastFunctionError]... % =compareBestChromosome(holdBestChromosome,holdLeastFunctionError,... % bestChromosome,leastFuntionError) %为染色体定义概率函数，好的染色体概率高，坏染色体概率低 %(6)p=chromosomeProbability(functionError); %按概率选择染色体函数： %(7)slecteChromosomeGroup=selecteChromome(fatherChromosomeGroup,p); %父代染色体杂交产生子代染色体函数 %(8)sonChrmosomeGroup=crossChromosome(slecteChromosomeGroup,2); %防止染色体超出解空间的函数 %(9)chromosomeGroup=checkSequence(chromosomeGroup,solutionSum) %变异函数 %(10)fatherChromosomeGroup=varianceCh(sonChromosomeGroup,0.8,solutionN); %通过实验有如下结果： %1。染色体应当多一些 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB 3-7变量非线性方程组的逆Broyden解法 matlab程序 function n_broyden clear all %清内存 clc %清屏 format long i=input('请输入非线性方程组个数(3-7)i='); switch i case 3 syms x y z x0=input('请输入初值（三维列向量[x;y;z])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(x,y,z)='); f2=input('请输入f2(x,y,z)='); f3=input('请输入f3(x,y,z)='); F=[f1;f2;f3]; J=jacobian(F,[x,y,z]); %求jacobi矩阵 Fx0=subs(F,{x,y,z},x0); Jx0=subs(J,{x,y,z},x0); H=inv(Jx0); x1=x0-H*Fx0; k=0; while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0; q=subs(F,{x,y,z},x1)-subs(F,{x,y,z},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1; x0=x1; Fx0=subs(F,{x,y,z},x0); x1=x1-H*Fx0; k=k+1; end x1 k case 4 syms a b c d x0=input('请输入初值（列向量[a;b;c;d])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(a,b,c,d)='); f2=input('请输入f2(a,b,c,d)='); f3=input('请输入f3(a,b,c,d)='); f4=input('请输入f4(a,b,c,d)='); F=[f1;f2;f3;f4]; J=jacobian(F,[a,b,c,d]); %求jacobi矩阵 Fx0=subs(F,{a,b,c,d},x0); Jx0=subs(J,{a,b,c,d},x0); H=inv(Jx0); x1=x0-H*Fx0; k=0; while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0; q=subs(F,{a,b,c,d},x1)-subs(F,{a,b,c,d},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1; x0=x1; Fx0=subs(F,{a,b,c,d},x0); x1=x1-H*Fx0; k=k+1; end x1 k case 5 syms a b c d e x0=input('请输入初值（列向量[a;b;c;d;e])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(a,b,c,d,e)='); f2=input('请输入f2(a,b,c,d,e)='); f3=input('请输入f3(a,b,c,d,e)='); f4=input('请输入f4(a,b,c,d,e)='); f5=input('请输入f5(a,b,c,d,e)='); F=[f1;f2;f3;f4;f5]; J=jacobian(F,[a,b,c,d,e]); %求jacobi矩阵 Fx0=subs(F,{a,b,c,d,e},x0); Jx0=subs(J,{a,b,c,d,e},x0); H=inv(Jx0); x1=x0-H*Fx0; 第7章 求解非线性方程 7.1 多项式运算在MATLAB 中的实现 一、多项式的表达 n 次多项式表达为：n a +??++=x a x a x a p(x )1-n 1-n 1n 0，是n+1项之和 在MATLAB 中，n 次多项式可以用n 次多项式系数构成的长度为n+1的行向量表示 [a0, a1,……an-1,an] 二、多项式的加减运算 设 有 两 个 多 项 式 n a +??++=x a x a x a p1(x )1-n 1-n 1n 0和 m b +??++=x b x b x b p2(x )1-m 1-m 1m 0。它们的加减运算实际上就是它们的对应系 数的加减运算。当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。 例2 计算()()1635223-+++-x x x x a=[1, -2, 5, 3]; b=[0, 0, 6, -1]; c=a+b 例3 设()6572532345++-+-=x x x x x x f ，()3532-+=x x x g ，求f(x)+g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; g1=[0, 0, 0, g];%为了和f 的次数找齐 f+g1, f-g1 三、多项式的乘法运算 conv(p1,p2) 例4 在上例中，求f(x)*g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; conv(f, g) 四、多项式的除法运算 [Q, r]=deconv(p1, p2) 表示p1除以p2，给出商式Q(x)，余式r(x)。Q,和r 仍为多项式系数向量 例4 在上例中，求f(x)/g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; [Q, r]=deconv(f, g) 五、多项式的导函数 p=polyder(P)：求多项式P 的导函数 p=polyder(P,Q)：求P·Q 的导函数 实验1.1 用matlab 求解线性方程组 第一节 线性方程组的求解 一、齐次方程组的求解 rref （A ） %将矩阵A 化为阶梯形的最简式 null （A ） %求满足AX ＝0的解空间的一组基，即齐次线性方程组的基 础解系 【例１】 求下列齐次线性方程组的一个基础解系，并写出通解： 我们可以通过两种方法来解： 解法1： >> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果： ans= 1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵，得化简后的方程 ??? ??=+--=+--=-+-0 22004321 43214321x x x x x x x x x x x x 取x2，x4为自由未知量，扩充方程组为 即 提取自由未知量系数形成的列向量为基础解系，记 所以齐次方程组的通解为 解法2： clear A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; B=null(A, 'r') % help null 看看加个‘r’是什么作用, 若去掉r ，是什么结果？ 执行后可得结果： B= 1 0 1 0 0 1 0 1 ?? ?=-=-0 04321x x x x ?????? ?====4 4432221x x x x x x x x ??? ??? ??????+????????????=????? ???????1100001142 4321x x x x x x , 00111????? ? ??????=ε, 11002????? ???????=ε2 211εεk k x += 中南大学 MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券 的自行百度。所需m文件照本文档做即可，即新建（FILE）→脚本（NEW-Sscript）→复制本文档代码→运行（会跳出 保存界面，文件名默认不要修改，保存）→结果。第 一题需要把数据文本文档和m文件放在一起。全部测 试无误，放心使用。本文档针对做牛顿法求非线性函 数题目的同学，当然第一题都一样，所有人都可以用。 ←记得删掉这段话 班级： 学号： 姓名： 一、《MATLAB程序设计实践》Matlab基础 表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散空间函数值来表示取向分布函数，是三维取向分布函数的一个实例。由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一个matlab程序画出如下的几种图形来分析其取向分布特征：（1）用Slice函数给出其整体分布特征； （2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25, 30, 35 … 90)切面上f分布情况(需要用到subplot函数)； (3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。 备注：数据格式说明 解： （1）将文件内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下： fid=fopen(''); for i=1:18 tline=fgetl(fid); end phi1=1;phi=1;phi2=1;line=0; f=zeros(19,19,19); while ~feof(fid) tline=fgetl(fid); data=str2num(tline); line=line+1; if mod(line,20)==1 phi2=(data/5)+1; phi=1; 数据说明部分，与作图无关 此方向表示f 随着φ1从0,5,10,15, 20 …到90的变化而变化 此方向表示f 随着φ从0,5,10,15, 20 …到90的变化而变化 表示以下数据为φ2＝0的数据，即f （φ1，φ，0） Matlab线性方程组求解 1. Gauss消元法： function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k); %计算行列式 end det=det*a(n,n); for k=n:-1:1 %回代求解 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(k)=b(k)/a(k,k); end Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 2. 列主元Gauss消去法： function x=detGauss(a,b) % Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0; %选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n 基于Matlab 实现牛顿迭代法解非线性方程组 已知非线性方程组如下 2211221212 10801080x x x x x x x ?-++=??+-+=?? 给定初值0(0,0)T x =，要求求解精度达到0.00001 首先建立函数F(x)，方程组编程如下，将F.m 保存到工作路径中： function f=F(x) f(1)=x(1)^2-10*x(1)+x(2)^2+8; f(2)=x(1)*x(2)^2+x(1)-10*x(2)+8; f=[f(1) f(2)]; 建立函数DF(x),用于求方程组的Jacobi 矩阵，将DF.m 保存到工作路径中： function df=DF(x) df=[2*x(1)-10,2*x(2);x(2)^2+1,2*x(1)*x(2)-10]; 编程牛顿迭代法解非线性方程组，将newton.m 保存到工作路径中： clear; clc x=[0,0]'; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',0,x(1),x(2)); N=4; for i=1:N y=df\f'; x=x-y; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',i,x(1),x(2)); if norm(y)<0.0000001 break ; else end end 运行结果如下： 0 0.0000000 0.0000000 1 0.8000000 0.8800000 2 0.9917872 0.9917117 3 0.9999752 0.9999685 4 1.0000000 1.0000000 本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（M üller ）法和迭代法的加速等）及其MATLAB 程序，求解非线性方程组的方法及其MATLAB 程序. 2.1 方程（组）的根及其MATLAB 命令 2.1.2 求解方程（组）的solve 命令 求方程f (x )=q (x )的根可以用MATLAB 命令: >> x=solve('方程f(x)=q(x)',’待求符号变量x ’) 求方程组f i (x 1,…,x n )=q i (x 1,…,x n ) (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 求解多项式方程（组）的roots 命令 如果)(x f 为多项式，则可分别用如下命令求方程0)(=x f 的根，或求导数)('x f （见表 2-1）. 2.1.4 求解方程（组）的fsolve 命令 如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots 命令.如果非线性方程（组）是含有超越函数，则无法使用roots 命令，需要调用MATLAB 系统中提供的另一个程序fsolve 来求解.当然，程序fsolve 也可以用于多项式方程（组），但是它的计算量明显比roots 命令的大. fsolve 命令使用最小二乘法（least squares method ）解非线性方程（组） (F X =)0 的数值解，其中X 和F （X ）可以是向量或矩阵.此种方法需要尝试着输入解X 的初始值（向量或矩阵）X 0，即使程序中的迭代序列收敛，也不一定收敛到(F X =)0的根（见例2.1.8）. fsolve 的调用格式： X=fsolve(F,X0) 输入函数)(x F 的M 文件名和解X 的初始值（向量或矩阵）X 0，尝试着解方程（组） 实用数值方法(Matlab) 综述报告题目：非线性方程（组）求根问题 小组成员 许多数学和物理问题归结为解函数方程f(x)=0。方程f(x)=0的解称为方程的根。对于非线性方程，在某个范围内往往不止一个根，而且根的分布情况可能很复杂，面对这种情况，通常先将考察的范围花费为若干个子段，然后判断哪些子段内有根，然后再在有根子段内找出满足精度要求的近似根。为此适当选取有根子段内某一点作为根的初始值近似，然后运用迭代方法使之足部精确化。这就是方程求根的迭代法。下面介绍书上的几种方法： 1、二分法 （1）方法概要： 假定函数f(x)在[a,b]上连续，且f(a)f(b)=0,则方程f(x)=0在[a,b]内一定有实根。取其中 将其二分，判断所求的根在的左侧还是右侧，得到一个新的有根区间 点 [],长度为[a,b]的一半。对新的有根区间继续实行上述二分手段，直至二分k次后有根区间[]长度 可见，如果二分过程无限继续下去，这些有限根区间最终必收敛于一点,该点就是所求的根。在实际计算过程中不可能完成这个无限过程，允许有一定的误差，则二分k+1次后 只要有根区间[]的长度小于，那么结果关于允许误差就能“准确”地满足方程f(x)=0。 (2)计算框图： 2、开方法 对于给定,求开方值 为此，可以运用校正技术设计从预报值生成校正值的迭代公式。自然希望校正值 能更好满足所给方程： 这是个关于校正量的近似关系式，如果从中删去二次项，即可化归为一次方程 解之有 从而关于校正值有如下开方公式 上述演绎过程表明，开方法的设计思想是逐步线性化，即将二次方程的求解画归为一次方程求解过程的重复。开方公式规定了预报值与校正值之间的一种函数关系,这里 为开方法的迭代函数。 3、Newton法 （1）方法概要 线性方程组求解 1.直接法 Gauss消元法： function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k); end det=det*a(n,n); for k=n:-1:1 %回代 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(k)=b(k)/a(k,k); end Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 列主元Gauss消去法： function x=detGauss(a,b) % Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0;% 选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n z=a(k,j);a(k,j)=a(r,j);a(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det; end 非线性方程求解 摘要：利用matlab软件编写程序，分别采用二分法、牛顿法和割线法求解非线性方程， 0 2= -x e x 的根，要求精确到三位有效数字，其中对于二分法，根据首次迭代结果，事先估计迭代次数，比较实际迭代次数与估计值是否吻合。并将求出的迭代序列用表格表示。对于牛顿法和割线法，至少取3组不同的初值，比较各自迭代次数。将每次迭代计算值求出，并列于表中。 关键词：matlab、二分法、牛顿法、割线法。 引言： 现实数学物理问题中，很多可以看成是解方程的问题，即f（x）=0的问题，但是除了极少简单方程的根可以简单解析出来。大多数能表示成解析式的，大多数不便于计算，所以就涉及到算法的问题，算法里面，具体求根时，一般先寻求根的某一个初始近似值，然后再将初始近似值逐步加工成满足精度要求为止，但是，我们知道，人为计算大大的加重了我们的工作量，所以大多用计算机编程，这里有很多可以计算的软件，例如matlab等等。 正文： 一、二分法 1 二分法原理：对于在区间[，]上连续不断且满足·＜0的函数， 通过不断地把函数的零点所在的区间一分为二，使区间的两个端点逐步逼近零点，进而得到零点近似值的方法叫做二分法。 2 二分法求根步骤：(1)确定区间，，验证·＜0，给定精确度；（2)求区间，的中点；(3)计算。若=，则就是函数的零点；若· ＜0，则令=；若·＜0，则令=。(4)判断是否达到精确度；即若 ＜，则得到零点近似值（或）；否则重复步骤2-4． 3 二分法具体内容：精度要求为5e-6，,解得实际迭代次数与估计值基本吻合，迭代如下表。n=2 c=0.000000 fc=-1.000000 n=11 c=-0.705078 fc=0.003065 n=3 c=-0.500000 fc=-0.356531 n=12 c=-0.704102 fc=0.001206 n= 4 c=-0.750000 fc=0.090133 n=13 c=-0.703613 fc=0.000277 n= 5 c=-0.625000 fc=-0.14463 6 n=14 c=-0.703369 fc=-0.00018 7 n=6 c=-0.687500 fc=-0.030175 n=15 c=-0.703491 fc=0.000045 n=7 c=-0.718750 fc=0.029240 n=16 c=-0.703430 fc=-0.000071 n= 8 c=-0.703125 fc=-0.000651 n=17 c=-0.703461 fc=-0.000013 n= 9 c=-0.710938 fc=0.014249 n=18 c=-0.703476 fc=0.000016MATLAB解线性方程组的直接方法
Matlab求解线性方程组非线性方程组
matlab程序设计实践-牛顿法解非线性方程
利用MATLAB求线性方程组
遗传算法解非线性方程组的Matlab程序
线性方程组求解matlab实现
3-7变量非线性方程组的逆Broyden解法matlab程序
MATLAB应用 求解非线性方程
实验一用matlab求解线性方程组
matlab程序设计实践-牛顿法解非线性方程
Matlab线性方程组求解(Gauss消去法)
基于Matlab的牛顿迭代法解非线性方程组
第二章非线性方程(组)的数值解法的matlab程序
MATLAB 非线性方程(组)求根
线性方程组求解Matlab程序(精.选)
用matlab对非线性方程求解
相关文档
- 用matlab对非线性方程求解
- 非线性方程组求解及matlab实现分解
- 用matlab求解非线性方程组的几种方法之程序
- MATLAB应用 求解非线性方程
- matlab程序设计实践-牛顿法解非线性方程
- 牛顿法求解非线性方程组matlab源程序
- Matlab求解线性方程组、非线性方程组
- 实验5 非线性方程求根及其MATLAB实现
- Matlab求解线性方程组非线性方程组
- MATLAB 非线性方程(组)求根
- matlab非线性方程的解法(含牛拉解法)
- MATLAB求解非线性方程
- 推荐-Broyden方法求解非线性方程组的Matlab实现 精品
- MATLAB应用 求解非线性方程
- 非线性方程组求解及matlab实现讲解
- matlab非线性方程求解
- 非线性方程组求解及matlab实现(精选)
- matlab程序设计实践-牛顿法解非线性方程
- 基于Matlab的牛顿迭代法解非线性方程组
- 实验2利用matlab解非线性、微分方程组答案