非线性方程组
非线性方程组的求解(汇编)
非线性方程组的求解摘要：非线性方程组求解是数学教学中，数值分析课程的一个重要组成部分，作为一门学科，其研究对象是非线性方程组。
求解非线性方程组主要有两种方法：一种是传统的数学方法，如牛顿法、梯度法、共轭方向法、混沌法、BFGS 法、单纯形法等。
传统数值方法的优点是计算精度高，缺点是对初始迭代值具有敏感性，同时传统数值方法还会遇到计算函数的导数和矩阵求逆的问题，对于某些导数不存在或是导数难求的方程，传统数值方法具有一定局限性。
另一种方法是进化算法，如遗传算法、粒子群算法、人工鱼群算法、差分进化算法等。
进化算法的优点是对函数本身没有要求，不需求导，计算速度快，但是精度不高。
关键字：非线性方程组、牛顿法、BFGS 法、记忆梯度法、Memetic 算法1： 三种牛顿法：Newton 法、简化Newton 法、修改的Newton 法【1-3】求解非线性方程组的Newton 法是一个最基本而且十分重要的方法, 目前使用的很多有效的迭代法都是以Newton 法为基础, 或由它派生而来。
n 个变量n 个方程的非线性方程组, 其一般形式如下:⎪⎪⎩⎪⎪⎨⎧===0),...,(...0),...,(0),...,(21212211n n n n x x x f x x x f x x x f （1）式(1)中,),...,(21n i x x x f ( i=1, ⋯, n) 是定义在n 维Euclid 空间Rn 中开域 D 上 的实值函数。
若用向量记号，令：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n x x x ...X 21,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡====)(...)()(0),...,(...0),..,(0)...,()(2121212,211X f X f X f x x x f x x x f x x x f X F nn n n n则方程组（1）也可表示为：0)(=X F（2） 其中：X ∈R n ,F ∶R n →R 0, F(X) ∈R n , R n 为赋值空间。
线性和非线性方程组的解法
3.2.4 迭代的收敛性
♦ 松弛法的收敛性分析； – A = Q -R = (D -ωL)/ω –[(1-ω)D+ωU]/ω； – QX = RX + B； – Xk+1 = Q-1RXk + Q-1B； – Xk+1 = (D -ωL)-1[(1-ω)D+ωU]Xk + Q-1B； ♦ 收敛的必要条件； – 0<ω<2；
3.2.2 Gauss-Seidel迭代法
♦ Gauss-Seidel迭代法；
k x1 +1 k+1 x2 k xn +1
k x1 +1 k+1 x2 k xn +1
= = L =
= =
[b −(a x +L+ a x )]/ a [b −(a x + a x +L+ a
Hale Waihona Puke [b −(a x +L+ a x )]/ a [b −(a x + a x +L+ a
1 k 12 2 k 21 1 k 1n n
11
]
X [k] = D−1B − (I − D−1A) X [k−1]
[
]
3.2.1 简单迭代法
♦ 简单迭代法的特征； – X[k] = D-1B - (I - D-1A)X[k-1] – 一阶；
L =
– Xk+1 = ω·(D-1LXk+1 + D-1UXk+D-1B)+(1-ω)·Xk； – Xk+1 = ω·(D-1LXk+1 + D-1UXk+D-1B)+(1-ω)·Xk；
第7章非线性方程组的数值解法
f 1 y f 2 2 y
2 y ( 1,1 ) 2
( 1,1 )
( y 3) ( 1, 1 )
( 1, 1 )
( x 1) ( 1 , 1 ) 2
( 1,1 )
f 1 f 2 2 2[ 2 * ( 3) ( 2 ) * ( 2 )] 4 f1 f2 g10 x ( 1,1) x ( 1,1) x f 1 f 2 g 2 2[ 2 * ( 3) 2 * ( 2 )] 20 20 y y f 1 y f 2 ( 1, 1 ) ( 1, 1 )
完
f ( x0 h, y0 k ) f ( x0 , y0 ) ( h k ) f ( x0 , y0 ) x y 1 2 ( h k ) f ( x 0 , y0 ) 2! x y 1 n ( h k ) f ( x 0 , y0 ) n! x y 1 n 1 ( h k ) f ( x0 h, y0 k ) ( n 1)! x y
2
2
令
0
得 f 1 f 1 ( g10 x g 20 y ) f 1 ( g10 ( g f 1 g f 1 ) 2 ( g 10 20 10 x y f 2 g 20 x f 2 g 20 x f 2 ) f2 y f 2 2 ) ( x y
1
f 1 ( x 0 , y0 ) f ( x , y ) 2 0 0
从n到n+1的迭代格式为：
f 1 ( x n , y n ) xn 1 x n x y y f 2 ( xn , yn ) n 1 n x
非线性代数方程组的解法
δ 1 = δ 0 + Δδ 1
11
(K1 )−1
=
Δδ 1 F1 − F0
= δ1 −δ0 ψ 1 −ψ 0
Δδ 2 = (K1)−1(R − F1)
………
(K i )−1 =
Δδ i Δψ i
=
δ ψ
i i
− δ i−1 −ψ i−1
(2.13)
显然 K i 就是相应于 Δδ i = δ i − δ i−1 与 Δψ i =ψ i −ψ i−1 的割线劲度。但实际上对于多维情
是对称矩阵。所以式（2.21）是对称秩 1 算法。 (2) 秩 2 算法
一个 N × N 阶的秩 2 矩阵，总可以表示为
[ ] 1T B2T
⎤ ⎥ ⎦
=
A1 B1T
+
A2 B2T
(2.22)
式中A1、A2、B1和B2均为N×1 维向量。将上式代入（2.14），再代入（2.15）得 A1B1T Δψi + A2 B2T Δψi = Δδi − (K i−1 )−1 Δψi
ψ(δi ) ≡ K (δi )δi − R ≠ 0
ψ(δ) 作为对平衡偏离的一种度量，称为失衡力。
对于一个单变量问题的非线性方程，直接迭代法的计算过程如图 2.1 和图 2.2 所示，它
们分别给出 F~δ为凸和凹曲线时的迭代过程。可以看出 K(δ)就是过曲线上点(δ, F(δ)与原点
的割线斜率。对于单变量问题，这一迭代过程是收敛的，但对多自由度情况，由于未知量通
似解。若
ψi−1 = ψ(δi−1) ≡ F (δi−1) − R ≠ 0
希望能找到一个更好的、方程（2.4）的近似解为
δ = δi = δi−1 + Δδi
解非线性方程组的牛顿迭代法
为克服这两个缺点，通常可用下述方法.
（1） 简化牛顿法，也称平行弦法.
xk 1 xk Cf ( xk )
其迭代公式为 （4.7）
C 0,1 ,.
迭代函数 ( x) x Cf ( x).
若在根 x * 附近成立 ( x) 1 Cf ( x) 1 ，即取 0 Cf ( x) 2，则迭代法（4.7）局部收敛.
8
xk
C 2 C
q2
k
1 q
2k
.
对任意 x0 0，总有 q 1，故由上式推知，当 k 时 xk C ，即迭代过程恒收敛. 例8 解 求 115 .
表7 6 计算结果 k 0 1 2 3 4 xk 10 10.750000 10.723837 10.723805 10.723805
f ( x) , f ( x)
由于
( x)
f ( x) f ( x) . 2 [ f ( x)]
假定 x *是 f ( x) 的一个单根，即 f ( x*) 0, f ( x*) 0 ， 则由上式知 ( x*) 0 ，于是依据定理4可以断定，牛顿法 在根 x *的邻近是平方收敛的.
准备 迭代
x0 ，计算 f 0 f ( x0 ), 选定初始近似值
步骤2
按公式
x1 x0 f 0 / f 0
迭代一次，得新的近似值 x1，计算 f1 f ( x1 ), f1 f ( x1 ). 步骤3 控制
x1 满足 1 如果
f1 2 ，则终 或
5
止迭代，以 x1作为所求的根；否则转步骤4. 允许误差，而
3
又因
( x*)
f ( x*) , f ( x*)
非线性方程组的数值方法
数值分析
二、简单迭代法
简单迭代法又称为不动点迭代法，基本思想是 首先构造不动点方程 x=φ(x)，即由方程 f(x)=0变换 为等价形式 x=φ(x), 式中φ(x)称为迭代函数。然后建 立迭代格式：xk+1 =φ(xk)称为不动点迭代格式。
当给定初值x0 后, 由迭代格式xk+1 =φ(xk)可求得数 列{xk}。如果{xk}收敛于α，且φ(x)在α连续，则α就是 不动点方程的根。因为：
lim
k
xk 1
lim
k
(
xk
)
(lim k
xk
)
知α=φ(α)，即{xk}收敛于方程的根 α 。
数值分析
数值分析
迭代法的几何意义
记y1=x , y2=φ(x) , 它们交点的横坐标α即为方程的根
y
y
y1 x
( x1, x1)
y2 (x)
( x2, x2)
( x0 ,( x0 ))
p
( x1,( x1 ))
③ 当k 时， xk 收敛到 x* ？ | x * xk | | ( x*) ( xk1) | | '(ξk1) | | x * xk1 |
L | x * xk1 | ...... Lk | x * x0 | 0
数值分析
数值分析
④
|
x * xk
|
L| 1 L
xk
xk1 |
xk+1
=
φ
(xk)
得到的序列
xk
k 0
收敛于φ(x) 在[a, b]上的唯一不动点。并且有误差估
计式：
|
x
*
xk
|
第六章 非线性方程(组)的求解
* * 又当 n 充 分 大[ 以 a ,b 后 ] ， (x ,x ), 于是 m 为偶数 n n 时， x [ a ,b ],f (x ) 0 ,不 变 号 了 ！ n n
2）二分法线性收敛，收敛因子为1/2。
* x x n 1 1 * 1 * x x ( x a ) ( x x ), . n n 1 n 1 n 1 * 2 2 x x2 n 1
f (x) m(x x*)m1h(x) h(x) g(x) 1 (x x*)g(x),h(x*) 0, m f (x) (x) x x 1 (x-x* )g(x) / h(x), f (x) m (x*) 1 1 , (x*) 1 1 1 , m m 牛顿迭代线性收敛，且 随 m的增加收敛性越来越差 。 重根时的改进：
或
定理一的条件太强，不便于实际应用。下面给出一个局部 收敛定理。
由迭代（ 6 -1 -1 ） 产产生的 x 均收 数 敛收敛 n * 1 x x x x n n n 1 1 L n L * 或 x x x x n 1 L1 0
* * * 定理二 :如果 (x ) 连续 (x ， ) L 1, 则 x N (x ， 0 δ )
关于初值的问题： 一般来说采用试探法，但对于一些实际问题初值的选择并 不困难，它是明确的。
关于重根的问题：
* 设 x 是 f( x ) 的 m 重零 m 点 1) , 此 （时 * m * f( x ) ( x x ) g ( x ), g ( x ) 0 , 1 * m 1 * f ( x ) m ( x x ) [ g ( x ) ( x x ) g ( x )], m
称算法（6-1-1）为牛顿迭代法。 f (x) 证明：令 (x) x ,则 xN (x*), f (x) 0 x (x) f (x) (x) f *) 0 (x) 1 f ( x ), ( x ,牛顿迭代收敛 2 [ f (x)] () * * 又 xn1 x (xn) (x ) (xn x*)2; 2 xn1 x* () c,至少二阶收敛。 2 2 (x x*)
matlab求解非线性方程组
非线性方程组求解1.mulStablePoint用不动点迭代法求非线性方程组的一个根function [r,n]=mulStablePoint(F,x0,eps)%非线性方程组：f%初始解：a%解的精度：eps%求得的一组解：r%迭代步数：nif nargin==2eps=1.0e-6;endx0 = transpose(x0);n=1;tol=1;while tol>epsr= subs(F,findsym(F),x0); %迭代公式tol=norm(r-x0); %注意矩阵的误差求法，norm为矩阵的欧几里德范数n=n+1;x0=r;if(n>100000) %迭代步数控制disp('迭代步数太多，可能不收敛！');return;endend2.mulNewton用牛顿法法求非线性方程组的一个根function [r,n]=mulNewton(F,x0,eps)if nargin==2eps=1.0e-4;endx0 = transpose(x0);Fx = subs(F,findsym(F),x0);var = findsym(F);dF = Jacobian(F,var);dFx = subs(dF,findsym(dF),x0);r=x0-inv(dFx)*Fx;n=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);dFx = subs(dF,findsym(dF),x0);r=x0-inv(dFx)*Fx; %核心迭代公式tol=norm(r-x0);n=n+1;if(n>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endend3.mulDiscNewton用离散牛顿法法求非线性方程组的一个根function [r,m]=mulDiscNewton(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=transpose(x0)-inv(J)*fx;m=1;tol=1;while tol>epsxs=r;fx = subs(F,findsym(F),xs);J = zeros(n,n);for i=1:nx1 = xs;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=xs-inv(J)*fx; %核心迭代公式tol=norm(r-xs);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endendformat short;4.mulMix用牛顿-雅可比迭代法求非线性方程组的一个根function [r,m]=mulMix(F,x0,h,l,eps)if nargin==4eps=1.0e-4;endn = length(x0);J = zeros(n,n);Fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));C =D - J;inD = inv(D);H = inD*C;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = Hm*inD*Fx;r = transpose(x0)-dr; m=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));C =D - J;inD = inv(D);H = inD*C;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = Hm*inD*Fx;r = x0-dr; %核心迭代公式tol=norm(r-x0);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endend5.mulNewtonSOR用牛顿-SOR迭代法求非线性方程组的一个根function [r,m]=mulNewtonSOR(F,x0,w,h,l,eps)if nargin==5eps=1.0e-4;endn = length(x0);J = zeros(n,n);Fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));L = -tril(J-D);U = -triu(J-D);inD = inv(D-w*L);H = inD*(D - w*D+w*L);;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = w*Hm*inD*Fx;r = transpose(x0)-dr;m=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));L = -tril(J-D);U = -triu(J-D);inD = inv(D-w*L);H = inD*(D - w*D+w*L);;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = w*Hm*inD*Fx;r = x0-dr; %核心迭代公式tol=norm(r-x0);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endend6.mulDNewton用牛顿下山法求非线性方程组的一个根function [r,m]=mulDNewton(F,x0,eps)%非线性方程组：F%初始解：x0%解的精度：eps%求得的一组解：r%迭代步数：nif nargin==2eps=1.0e-4;endx0 = transpose(x0);dF = Jacobian(F);m=1;tol=1;while tol>epsttol=1;w=1;Fx = subs(F,findsym(F),x0);dFx = subs(dF,findsym(dF),x0);F1=norm(Fx);while ttol>=0 %下面的循环是选取下山因子w的过程r=x0-w*inv(dFx)*Fx; %核心的迭代公式Fr = subs(F,findsym(F),r);ttol=norm(Fr)-F1;w=w/2;endtol=norm(r-x0);m=m+1;x0=r;if(m>100000) %迭代步数控制disp('迭代步数太多，可能不收敛！');return;endend7.mulGXF1用两点割线法的第一种形式求非线性方程组的一个根function [r,m]=mulGXF1(F,x0,x1,eps)format long;if nargin==3eps=1.0e-4;endx0 = transpose(x0);x1 = transpose(x1);n = length(x0);fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);for i=1:nxt = x1;xt(i) = x0(i);J(:,i) = (subs(F,findsym(F),xt)-fx1)/h(i);endr=x1-inv(J)*fx1;m=1;tol=1;while tol>epsx0 = x1;x1 = r;fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);for i=1:nxt = x1;xt(i) = x0(i);J(:,i) = (subs(F,findsym(F),xt)-fx1)/h(i);endr=x1-inv(J)*fx1;tol=norm(r-x1);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endendformat short;8.mulGXF2用两点割线法的第二种形式求非线性方程组的一个根function [r,m]=mulGXF2(F,x0,x1,eps)format long;if nargin==3eps=1.0e-4;endx0 = transpose(x0);x1 = transpose(x1);n = length(x0);fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);xt = x1;xt(1) = x0(1);J(:,1) = (subs(F,findsym(F),xt)-subs(F,findsym(F),x1))/h(1);for i=2:nxt = x1;xt(1:i) = x0(1:i);xt_m = x1;xt_m(1:i-1) = x0(1:i-1);J(:,i) = (subs(F,findsym(F),xt)-subs(F,findsym(F),xt_m))/h(i);endr=x1-inv(J)*fx1;m=1;tol=1;while tol>epsx0 = x1;x1 = r;fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);xt = x1;xt(1) = x0(1);J(:,1) = (subs(F,findsym(F),xt)-subs(F,findsym(F),x1))/h(1);for i=2:nxt = x1;xt(1:i) = x0(1:i);xt_m = x1;xt_m(1:i-1) = x0(1:i-1);J(:,i) = (subs(F,findsym(F),xt)-subs(F,findsym(F),xt_m))/h(i);endr=x1-inv(J)*fx1;tol=norm(r-x1);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endendformat short;9.mulVNewton用拟牛顿法求非线性方程组的一组解function [r,m]=mulVNewton(F,x0,A,eps)%方程组：F%方程组的初始解：x0% 初始A矩阵：A%解的精度：eps%求得的一组解：r%迭代步数：mif nargin==2A=eye(length(x0)); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendx0 = transpose(x0);Fx = subs(F, findsym(F),x0);r=x0-A\Fx;m=1;tol=1;while tol>epsx0=r;Fx = subs(F, findsym(F),x0);r=x0-A\Fx;y=r-x0;Fr = subs(F, findsym(F),r);z= Fr-Fx;A1=A+(z-A*y)*transpose(y)/norm(y); %调整A A=A1;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endtol=norm(r-x0);end10.mulRank1用对称秩1算法求非线性方程组的一个根function [r,n]=mulRank1(F,x0,A,eps)if nargin==2l = length(x0);A=eye(l); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-inv(A)*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-inv(A)*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;A1=A+ fr *transpose(fr)/(transpose(fr)*y); %调整A A=A1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多，可能不收敛！');return;endtol=norm(r-x0);end11.mulDFP用D-F-P算法求非线性方程组的一组解function [r,n]=mulDFP(F,x0,A,eps)if nargin==2l = length(x0);B=eye(l); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-B*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-B*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;B1=B+ y*y'/(y'*z)-B*z*z'*B/(z'*B*z); %调整AB=B1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多，可能不收敛！');return;endtol=norm(r-x0);end12.mulBFS用B-F-S算法求非线性方程组的一个根function [r,n]=mulBFS(F,x0,B,eps)if nargin==2l = length(x0);B=eye(l); %B取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-B*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-B*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;u = 1 + z'*B*z/(y'*z);B1= B+ (u*y*y'-B*z*y'-y*z'*B)/(y'*z); %调整B B=B1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多，可能不收敛！');return;endtol=norm(r-x0);end13.mulNumYT用数值延拓法求非线性方程组的一组解function [r,m]=mulNumYT(F,x0,h,N,eps)format long;if nargin==4eps=1.0e-8;endn = length(x0);fx0 = subs(F,findsym(F),x0);x0 = transpose(x0);J = zeros(n,n);for k=0:N-1fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endinJ = inv(J);r=x0-inJ*(fx-(1-k/N)*fx0);x0 = r;endm=1;tol=1;while tol>epsxs=r;fx = subs(F,findsym(F),xs);J = zeros(n,n);for i=1:nx1 = xs;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=xs-inv(J)*fx; %核心迭代公式tol=norm(r-xs);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endendformat short;14.DiffParam1用参数微分法中的欧拉法求非线性方程组的一组解function r=DiffParam1(F,x0,h,N)%非线性方程组：f%初始解：x0%数值微分增量步大小：h%雅可比迭代参量:l%解的精度：eps%求得的一组解：r%迭代步数：nx0 = transpose(x0);n = length(x0);ht = 1/N;Fx0 = subs(F,findsym(F),x0);for k=1:NFx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endinJ = inv(J);r = x0 - ht*inJ*Fx0;x0 = r;end15.DiffParam2用参数微分法中的中点积分法求非线性方程组的一组解function r=DiffParam2(F,x0,h,N)%非线性方程组：f%初始解：x0%数值微分增量步大小：h%雅可比迭代参量:l%解的精度：eps%求得的一组解：r%迭代步数：nx0 = transpose(x0);n = length(x0);ht = 1/N;Fx0 = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nxt = x0;xt(i) = xt(i)+h(i);J(:,i) = (subs(F,findsym(F),xt)-Fx0)/h(i);endinJ = inv(J);x1 = x0 - ht*inJ*Fx0;for k=1:Nx2 = x1 + (x1-x0)/2;Fx2 = subs(F,findsym(F),x2);J = zeros(n,n);for i=1:nxt = x2;xt(i) = xt(i)+h(i);J(:,i) = (subs(F,findsym(F),xt)-Fx2)/h(i);endinJ = inv(J);r = x1 - ht*inJ*Fx0;x0 = x1;x1 = r;end16.mulFastDown用最速下降法求非线性方程组的一组解function [r,m]=mulFastDown(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endlamda = fx/sum(diag(transpose(J)*J));r=x0-J*lamda; %核心迭代公式fr = subs(F,findsym(F),r);tol=dot(fr,fr);x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endendformat short;17.mulGSND用高斯牛顿法求非线性方程组的一组解function [r,m]=mulGSND(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endDF = inv(transpose(J)*J)*transpose(J);r=x0-DF*fx; %核心迭代公式tol=norm(r-x0);x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endendformat short;18.mulConj用共轭梯度法求非线性方程组的一组解function [r,m]=mulConj(F,x0,h,eps)format long;if nargin==3eps=1.0e-6;endn = length(x0);x0 = transpose(x0);fx0 = subs(F,findsym(F),x0);p0 = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)*(1+h);p0(:,i) = -(subs(F,findsym(F),x1)-fx0)/h;endm=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endlamda = fx/sum(diag(transpose(J)*J));r=x0+p0*lamda; %核心迭代公式fr = subs(F,findsym(F),r);Jnext = zeros(n,n);for i=1:nx1 = r;x1(i) = x1(i)+h;Jnext(:,i) = (subs(F,findsym(F),x1)-fr)/h;endabs1 = transpose(Jnext)*Jnext;abs2 = transpose(J)*J;v = abs1/abs2;if (abs(det(v)) < 1)p1 = -Jnext+p0*v;elsep1 = -Jnext;endtol=norm(r-x0);p0 = p1;x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');return;endendformat short;19.mulDamp用阻尼最小二乘法求非线性方程组的一组解function [r,m]=mulDamp(F,x0,h,u,v,eps)format long;if nargin==5eps=1.0e-6;endFI = transpose(F)*F/2;n = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsj = 0;fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;afx = subs(F,findsym(F),x1);J(:,i) = (afx-fx)/h;endFIx = subs(FI,findsym(FI),x0);for i=1:nx2 = x0;x2(i) = x2(i)+h;gradFI(i,1) = (subs(FI,findsym(FI),x2)-FIx)/h;ends=0;while s==0A = transpose(J)*J+u*eye(n,n);p = -A\gradFI;r = x0 + p;FIr = subs(FI,findsym(FI),r);if FIr<FIxif j == 0u = u/v;j = 1;elses=1;endelseu = u*v;j = 1;if norm(r-x0)<epss=1;endendendx0 = r;tol = norm(p);m=m+1;if(m>100000) %迭代步数控制disp('迭代步数太多，可能不收敛！');return;endendformat short;。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
非线性方程组前言非线性方程组，顾名思义就是未知数的幂除了不是1，其他都有可能！线性方程组其实只是非常小的一类，非线性方程组才是大类！也正因此非线性方程组包含各种各样的方程形式，所以它的解和对应的求解方法不可能像线性方程组那样完美，即都是局部收敛的。
先给出一个直观的非线性方程组例子：个人对两个问题的理解：1、非线性方程组如果有解，一般都有很多解！如何理解：把方程组的解看成是各个函数图像的交点的。
我们知道非线性方程组的各个函数就都是复杂曲线、面，甚至是高纬空间里的复杂东西；线性方程组的各个函数就是最简单的直线、面！各个复杂函数图像间的相交机会很多，并且只要相交，就是多个交点(因为交线、交面里有无数的交点)，也就是有多个解，可以想象，非线性方程组有多解是很平常的一件事，对于复杂的非线性函数没解才不正常！可以想象，这些解是等价的！没有说是等级更高，谁等级低一些。
都是解！因为：只要是解，它就只满足一个条件：让方程组中的各个方程=0。
所以无法用什么评判标准(比如范数)来说哪个解的等级高一些或者效果更好一些。
注意：这里的解等价和欠定线性方程组通解中的唯一极小范数解不一样！可以想象二者的区别：非线性方程组中的解都是实打实存在的；而欠定线性方程组中除了特解，其他通解中的解说存在也行，说不存在那就是因为方程条件(个数)都不够！这些是啥都行的通解和非线性方程组中实打实存在的解肯定不能比！这样的话各个非线性方程组的局部收敛性就可以理解，即：空间中有很多解时，我每次只能找一个，那我找谁？找离我出发点最近的那个解呗。
所以不同的出发点，就有可能找到不同的解，这就是局部收敛性。
意思是：每个方程中把所有和有关的用一个变量代替，所有有关的用一个变量代替，即方程1中用：，但是很明显方程2的第一项两个变量相乘，没法用变量代替，并且，即使在方程2中能代替，那么就会有和，这样总未知数变成4个而方程只有2个，还是解不了。
所以，非线性方程组不可能用简单的线性变量代换来解。
本文介绍最常用、最好用的非线性方程组解法，包括牛顿法和拟牛顿法(割线法)两大类：牛顿法：原始牛顿法、修正牛顿法；拟牛顿法：逆Broyden秩1法、逆Broyden秩1第二公式、BFS秩2法；上面这5种方法的函数表达式、计算流程、代码都会详细说明和展示。
并且，这5种方法都很经济(算的快)、实用(适用各种非刁钻的函数)、易用(计算流程好理解，便于编程)。
本文侧重的是方法的使用，不提推导和太具体的其他细节。
本文要介绍的5种方法可分为两大类：牛顿法、拟牛顿法。
先把两类方法的优缺点列出来：牛顿法：用jacobi矩阵(导数)优点：导数法收敛速度巨快(平方收敛)；自校正误差不会传递；缺点：求导稍费事；只要赋值后的jacobi矩阵存在稀疏性、奇异性、病态等，就跪了；拟牛顿法：割线法思想，用近似矩阵趋近jacobi矩阵优点：jacobi矩阵的问题在这里都不是问题！这个优点极大提高解法的稳定性缺点：收敛速度介于平方收敛和直线收敛之间，稍慢一丢丢。
其实，牛顿法看上去要求导很麻烦，其实导数就求了一次而已！剩下都是在循环里对jacobi矩阵的赋值！所以去求导其实不是大问题。
大问题就是每次赋值后的jacobi矩阵要求逆！！这就对数值矩阵的要求很高了！并且实际问题中经常出现赋值后的jacobi矩阵是稀疏的。
所以，如果原函数们不刁钻，两种方法都可胜任。
但如果稍微感觉原函数有些复杂时，建议牺牲一丢的速度选择拟牛顿法！拟牛顿法的稳定性真的提高一个量级。
下面我们将正式开始方法的介绍，给定一个的非线性方程组的通式：对应的jacobi矩阵为：原始牛顿法自己给定，然后进行迭代：上述迭代可以用了，但是每次循环都要求逆很慢！所以换成下面这种形式：设真实解为，实现流程如下：步1：给出初始近似及计算精度和；步2：假定已进行了k次迭代，已求出和，计算并做赋值：步3：解上面的线性方程组(用)：步4：求下面两个式子：步5：若或，则置：并转步6；否则进行如下赋值后回到步2：步6：结束。
针对最开始的例子，下面给出matlab程序：clear ; clc% 未知数:ym 1 2;% 方程组: 统一用列向量f1 = 1^2 - 101 + 2^2 + 8;f2 = 12^2 + 1 - 102 + 8;% f1 = 41 - 2 + 0。
1ep(1)-1;% f2 = -1 + 42 + 1、81^2;= [1;2];f = [f1;f2];% 初始值: 统一用列向量0 = [2;3];error_dk = double( input('dk范数的精度:') );error_fkk = double( input('fkk范数的精度:') );num = input('停止迭代次数:');% jacobi1 = [diff(f1,1) diff(f1,2);diff(f2,1) diff(f2,2)] % 直接用自带函数求雅克比矩阵:jacobi = jacobian([f1,f2],[1,2]);for k = 1:numAk = double( ub(jacobi, , 0) );bk = double( ub(f, , 0) );dk = pre_eidel(Ak,-bk,k); % 步长0 = 0 + dk;fkk = double( ub(f, , 0) ); % fk+1单纯用来判断if norm(dk) < error_dk ， norm(fkk) < error_fkkbreak;endendif k < num_reult = 0;fprintf('精度已达要求，迭代提前结束!\n');fprintf('%d次迭代后, 近似解为:\n',k);_reultele_reult = 0;fprintf('迭代次数已达上限!\n');fprintf('似解为:\n',k);_reultendfprintf('f1结果为:%f\n',double( ub(f1,,0) ));fprintf('f2结果为:%f\n',double( ub(f2,,0) ));fprintf('dk范数:%f\n',norm(dk));fprintf('fkk范数:%f\n',norm(fkk));结果：dk范数的精度:0。
0001fkk范数的精度:0。
0001停止迭代次数:200第1次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0。
8306第2次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0。
0668第3次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0。
1969第4次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0。
2209第5次求解线性方程组, 当前万能赛德尔迭代矩阵谱半径为(越小越好): 0。
2214精度已达要求，迭代提前结束!5次迭代后, 近似解为:_reult =0。
99991、0000f1结果为:0。
000728f2结果为:0。
000184dk范数:0。
000000fkk范数:0。
000751修正牛顿法对于上面的原始牛顿法，如果每步计算改为固定的可得：这样就变成了"简化牛顿法"。
很明显可以看到简化牛顿法是线性收敛的！计算量小，但是收敛速度非常慢。
如果既拥有"原始牛顿法"的收敛速度，又拥有"简化牛顿法"的工作量省？，修正牛顿法。
对应的迭代格式为：从迭代公式可以得知，在每一个k的大循环内都有一个从1到m的小循环来求，下面同样给出实现流程：步1：给出初始近似及计算精度和；步2：根据方程阶数、个数n，求出效率最大化的m值或人为给定一个m值；步3：假定已进行k此迭代，已算出和，计算并赋值：：进入每次的内层循环，假设已内层循环i次，已算出和，先做1个赋值，再做2个计算：：先做3个赋值然后做一个计算：然后再一次判断：如果，转到步6(出小循环)，否则回到步4(没出小循环)；步6：若或，则赋值并结束：否则，然后重回步3(大循环)。
仍针对最开始的例子，给出matlab程序：clear ; clc;% 未知数:ym 1 2;% 方程组: 统一用列向量f1 = 1^2 - 101 + 2^2 + 8;f2 = 12^2 + 1 - 102 + 8;= [1;2];f = [f1;f2];% 初始值: 统一用列向量0 = [5;4];error_z = double( input('z范数的精度:') );error_fk = double( input('fk范数的精度:') );num = input('停止迭代次数:');% 直接用自带函数求雅克比矩阵:jacobi = jacobian([f1,f2],[1,2]);% 小循环m取多少的判断: 和方程个数N有关,求w的最大值ym M N;mn0 = [N;M];% 参数zn是方程的个数:zn = length();% 原始的效率对比方程:w = (N+1)log(M+1)、( (N+M)log(2) );% 让w方程求最大值的zm值:zm = double( olve( diff( ub(w,N,zn),M ) ) );mn1 = [zn;zm];% 最高效率值:wa = double( ub(w,mn0,mn1) );% 开始修正牛顿迭代:ki = 0;for k = 1:numfk = double( ub(f,,0) );Ak = inv( double( ub(jacobi, , 0) ) );for m = 1:round(zm)b = fk;0 = 0 - Akb;fki = double( ub(f,,0) );fk = fki;z = Akb;endif norm(z) < error_z ， norm(fk) < error_fk break;endendif k < num_reult = 0;fprintf('精度已达要求，迭代提前结束!\n'); fprintf('%d次迭代后, 近似解为:\n',k);_reultele_reult = 0;fprintf('迭代次数已达上限!\n');fprintf('似解为:\n',k);_reultendfprintf('f1结果为:%f\n',double( ub(f1,,0) ));fprintf('f2结果为:%f\n',double( ub(f2,,0) ));fprintf('z范数:%f\n',norm(z));fprintf('fk范数:%f\n',norm(fk));效果和原始牛顿法一样，可以感觉到速度提升了！！这里用到了求最效率的m值，下面对它的求法加以补充(完全可以不求手动给)：修正牛顿法与原始牛顿法的效率之比为：其中n就是方程组中方程的个数，这对于每个方程组都是固定常数。