不求导数矩阵的非线性方程组算法

1. matlab解非线性方程组

1. fsolve
求解非线性方程组
方程：
F(x)=0
x是一个向量，F(x)是该向量的函数向量，返回向量值
2. 语法
x = fsolve(fun,x0)
x = fsolve(fun,x0,options)
[x,fval] = fsolve(fun,x0)
[x,fval,exitflag] = fsolve(...)
[x,fval,exitflag,output] = fsolve(...)
[x,fval,exitflag,output,jacobian] = fsolve(...)
3. 描述
fsolve用于寻找非线性系统方程组的零点。
x = fsolve(fun,x0)以x0为初始值，努力寻找在fun中描述的方程组。
x = fsolve(fun,x0,options) 以x0为初始值，按照指定的优化设置“options”努力寻找在fun中描述的方程组。使用optimset设置这些选项。
[x,fval] = fsolve(fun,x0)返回在解x处的目标函数fun的值
[x,fval,exitflag] = fsolve(...)返回exitflag表示退出条件。
[x,fval,exitflag,output] = fsolve(...)返回output结构，该结构包含了优化信息。
[x,fval,exitflag,output,jacobian] = fsolve(...)返回在解x处的Jacobian函数。
4. 输入参数
4.1. fun
非线性系统方程。它是一个函数，以x作为输入，返回向量F。函数fun可以被指定为一个M文件函数的函数句柄。
x = fsolve(@myfun,x0)
这里的myfun是一个matlab函数，形如：
function F = myfun(x)
F = ... % Compute function values at x
fun也可以是一个异步函数的函数句柄：
x = fsolve(@(x)sin(x.*x),x0);
若用户定义的值为矩阵，则会被自动转换为向量。
若Jacobian能被计算出来且通过options = optimset('Jacobian','on')设置Jacobian选项为”on”，则函数fun必须在第2个输出参数中返回x处的Jacobian值J（它是一个矩阵）。注意：通过检查nargout的值，当fun被只带一个输出参数调用的情况下，该函数可避免计算J，仅只有一个输出参数。（这种情况下，优化算法仅需要知道F而不需J）。
function [F,J] = myfun(x)
F = ... % objective function values at x
if nargout > 1 % two output arguments
J = ... % Jacobian of the function evaluated at x
end
4.2. options
提供该函数有关的特定信息。
5. 输出参数
5.1. exitflag
一个用来表示算法终止原因的整数。
1：函数收敛到x
2：x的变化已经处在容许范围内
3：残差变化已经处在容许范围内
4：搜索方向飞幅度比指定的误差小
5：迭代次数超过options.MaxIter或函数估值的次数超过options.FunEvals
-1：算法被输出函数终止
-2：算法似乎收敛到一个非根点。
-3：可信半径变得太小
-4：沿当前方向的线性搜索不能足够地减小残差
5.2. output
包含关于优化信息的一个结构，其具有如下字段：
iterations：已经迭代的次数
funcCount：函数估值的次数
algorithm：所使用的算法
cgiterations：PCG迭代次数（仅适用于大规模算法）
stepsize：最终采取的步长（仅适用于中等规模算法）
firstorderopt：第1阶优化的观测值 。
5.3. options
优化设置。一些选项设置用于所有算法，部分与大规模算法（large-scale algrithm）。相关，部分与中等规模算法相关。可以使用optimset改变其中的设置。LargeScale选项指定使用哪种算法。
设为‘on’使用大规模算法，设为‘off’使用中等规模算法。
5.3.1. Medium-Scale and Large-Scale Algorithms
如下选项用于大规模和中等规模算法：
DerivativeCheck：将用户提供导数与有限差分导数相比较
Diagnostics：显示被解函数的诊断信息
DiffMaxChange：有限差分变量中的最大变化
DiffMinChange：有限差分变量中的最小变化
Display：显示的级别，‘off‘不显示输出，’iter‘显示每一步迭代的输出，’final‘显示最终的输出（默认）
FunValCheck：检查目标函数值是否有效。设为‘on’，当函数返回值为复数、inf或NaN将返回一个错误，设为‘off’将不显示错误。
Jacobian：设为‘ON’，fsolve将使用用户定义的Jacobian或Jacobian信息来估值目标函数，若设为‘off’，则使用有限差分逼近Jacobian。
MaxFunEvals：最大允许估值次数
MaxIter：最大迭代次数
OutputFcn：指定一个或多个输出函数，优化函数在每一个迭代过程中将调用这些函数。
PlotFcns：算法执行时显示进度条。从预定义中选择或自定义进度条。指定@optimplotx显示当前的点，@optimplotfunccount打印出函数的计数，@optimplotfval打印出函数值，@optimplotresnorm打印出残差范数，@optimplotstepsize打印出步长，@optimplotfirstorderopt打印优化参数的第1阶。
TolFun：函数值的终止误差。
TolX：x的终止误差
5.3.2. 仅适用于Large-Scale Algorithm
JacobMult：Jacobian乘法函数的函数句柄
JacobPattern
MaxPCGIter
PrecondBandWidth
TolPCG
5.3.3. 仅适用于Medium-Scale Algorithm
NonlEqnAlgorithm：
'dogleg' — Trust-region dogleg algorithm
(default)
'lm' — Levenberg-Marquardt
'gn' — Gauss-Newton
LineSearchType：'lm' (Levenberg-Marquardt)
'gn' (Gauss-Netwton) algorithms.

6. 使用优化工具箱完成以上函数操作
命令：optimtool

2. 非线性方程组的解法matlab

用matlab求解非线性方程组方法，可以用下列方法来实现：
方法一，使用solve函数求解
x = optimvar('x');
y = optimvar('y');
prob = optimproblem;
prob.Objective = -x - y/3;
prob.Constraints.cons1 = x + y <= 2;
prob.Constraints.cons2 = x + y/4 <= 1;
prob.Constraints.cons3 = x - y <= 2;
prob.Constraints.cons4 = x/4 + y >= -1;
prob.Constraints.cons5 = x + y >= 1;
prob.Constraints.cons6 = -x + y <= 2;
sol = solve(prob)

方法二，使用fsolve函数求解
F = @(x) [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))];

x0=[-5;-5];

[x,fval] = fsolve(F,x0)

方法三，使用迭代法求解，如Newton迭代法
m=3;
x0=zeros(m,1);
tol=1e-6;
x=x0-dfun(x0)\fun(x0);
n=1;
while(norm(x-x0>tol)) & n<1000
x0=x;
x=x0-dfun(x0)\fun(x0);
n=n+1;
end
x
这里，fun是原方程组，dfun是导数方程组

3. 如何解矩阵非线性方程组

只有1个矩阵是不能转化成非线性方程组的，还要知道一个列向量，这个列向量作为非其次方程组等号右边的量
然后矩阵的每行的每个元素后面依次加上x1...xn(如果这行有n个元素的话，也就是方程组有n个变量)，第2行也是这样，一直到第n行（假设有n个方程组）

4. 非线性方程组数值解法的介绍

20世纪60年代中期以后，发展了两种求解非线性方程组(1)的新方法。一种称为区间迭代法或称区间牛顿法，它用区间变量代替点变量进行区间迭代，每迭代一步都可判断在所给区间解的存在惟一性或者是无解。这是区间迭代法的主要优点，其缺点是计算量大。另一种方法称为不动点算法或称单纯形法，它对求解域进行单纯形剖分，对剖分的顶点给一种恰当标号，并用一种有规则的搜索方法找到全标号单纯形,从而得到方程(1)的近似解。这种方法优点是，不要求f(□)的导数存在,也不用求逆，且具有大范围收敛性，缺点是计算量大

5. 非线性方程组数值解法的割线法

若对方程组 (1)线性化时使用插值方法确定线性方程组

(6)
中的Ak和bk,则可得到一类称为割线法的迭代序列。假定已知第k步近似尣k，为确定Ak和bk,可在尣附近取n个辅助点у忋(j=1,2,…,n),使n个向量

线性无关，由插值条件可知

由此可求得

由(6)解得

以此作为方程 (1)的新近似，记作

，于是得到

(7)
(7)称为解非线性方程组的割线法。辅助点у忋 取得不同就得到不同的割线法程序，例如取

为常数(j=1,2,…,n)，就得到与(5)相同的程序，由于它只依赖于尣点的信息，故也称一点割线法,若取

它依赖于点尣及

, 称为两点割线法。其他多点割线法由于稳定性差，使用较少。
非线性方程组数值解法 - 布朗方法 布朗采用对每个分量方程 ƒi(尣)=0逐个进行线性化并逐个消元的步骤，即在每迭代步中用三角分解求线性方程组的解，得到了一个效率比牛顿法提高近一倍的迭代法，即

式中

(8)中当i=n时求得xn记作

,再逐次回代,求出

(i=n-1,n-2,…,1)就完成了一个迭代步。布朗迭代程序的敛速仍保持p=2，而每一迭代步的工作量

，故效率

对这方法还可与牛顿法一样进行改进，得到一些效率更高的算法。这类方法是70年代以来数值软件包中常用的求解非线性方程组的算法。
非线性方程组数值解法 - 拟牛顿法 为减少牛顿法的计算量，避免计算雅可比矩阵及其逆，60年代中期出现了一类称为拟牛顿法的新算法,它有不同的形式，常用的一类是秩1的拟牛顿法，其中不求逆的程序为

式中

,

,

，称为逆拟牛顿公式。计算时先给出尣及 B0，由(9)逐步迭代到满足精度要求为止。每步只算 n个分量函数值及O(n)的计算量,比牛顿法一步计算量少得多。理论上已证明，当尣及B0选得合适时,它具有超线性收敛速度,但实践表明效率并不高于牛顿法，理论上尚无严格证明。
非线性方程组数值解法 - 最优化方法 求方程组 (1)的问题等价于求目标函数为

的极小问题,因此可用无约束最优化方法求问题(1)的解（见无约束优化方法）。 非线性方程组数值解法 - 连续法 又称嵌入法，它可以从任意初值出发求得方程组(1)的一个足够好的近似解,是一种求出好的迭代初值的方法。连续法的基本思想是引入参数 t∈【0,b】，构造算子H(尣,t)，使它满足条件：H(尣,0)=ƒ0(尣),H(尣,b)=ƒ(尣)，其中ƒ0(尣)=0的解尣0是已知，方程：

(10)
在t∈【0,b】上有解尣=尣(t)，则尣(b)=尣就是方程(1)的解。当b有限时，通常取b=1，例如可构造。

(11)
这里尣是任意初值,显然H(尣,0)=0,H(尣,1)=ƒ(尣)。为了求得(10)在t=1的解尣=尣(1),可取分点0=t0<t1<…<tN=1在每个分点ti(i=1,2,…,N)上，求方程组 H(尣,ti)=0 (i=1,2,…,N) (12)
的解尣，如果取尣为初值，只要

足够小，牛顿迭代就收敛，但这样做工作量较大。已经证明，如果方程组(12)只用一步牛顿法，当t=tN=1时，再用牛顿迭代，结果仍具有2阶收敛速度。以(11)为例,得到连续法的程序为：

若H(尣,t)的偏导数Ht(尣,t)及

在D×【0,1】嶅R

上连续。且

非奇异，则由(10)对t求导可得到等价的微分方程初值问题:

(13)
于是求方程(10)的解就等价于求常微初值问题(13)的解，求(13)的解可用数值方法由t=0计算到t=tN=b得到数值解

。已经证明只要N足够大,以尣为初值再进行牛顿迭代可收敛到方程(1)的解x，这种算法称为参数微分法。
20世纪60年代中期以后，发展了两种求解非线性方程组(1)的新方法。一种称为区间迭代法或称区间牛顿法，它用区间变量代替点变量进行区间迭代，每迭代一步都可判断在所给区间解的存在惟一性或者是无解。这是区间迭代法的主要优点，其缺点是计算量大。另一种方法称为不动点算法或称单纯形法，它对求解域进行单纯形剖分，对剖分的顶点给一种恰当标号，并用一种有规则的搜索方法找到全标号单纯形,从而得到方程(1)的近似解。这种方法优点是，不要求ƒ(尣)的导数存在,也不用求逆，且具有大范围收敛性，缺点是计算量大。

6. 非线性方程组解法

非线性方程组可以表示为：

在位移为基本未知量的有限元分析中， 是结点位移向量， 是结点载荷向量。对于线性方程组 ，由于 是常数矩阵，可以没有困难直接求解。对于非线性方程组需要迭代求解。

只适用于与变形历史无关的非线性问题，例如非线性弹性问题，利用形变理论分析的弹塑性问题，稳态蠕变问题等。对于依赖于变形历史的非线性问题，应力需要由应变所经历的路径决定，直接迭代法不适用。例如加载路径不断变化或涉及卸载和反复加载等弹塑性问题。此时要利用增量理论。
可以指出，当 是凸曲线（如图所示），其中 是标量，即系统为单自由度情形，通常解收敛。

(1)方程可以改写为
选择一个初始的试探解 代入上式可以得到初始的 ，接着可以得到第一步迭代的位移解 ，由此类推得到第 次迭代后的近似解 ，一直到误差的某种范数小于某个规定的容许小量 ，即 ，上述迭代终止。

为避免每次迭代对刚度矩阵 进行求逆计算，一般可以将刚度矩阵设定为常数进行迭代过程，单自由度系统下的常刚度直接迭代法如图所示

一般情况下，（1）不能被精确地满足，即 ，为了得到进一步的近似解 （假设 为某单自由度方向上的位移） ，将 表示成在 附近的仅保留线性项的Taylor展开式

且有
式中 是切线矩阵(tangent matrix)，即



其迭代过程如图所示

核心思想是将载荷分解成若干增量步，即 ，相应的位移也分为同样的步数，即 ，每两步之间的增长量称之为增量。增量解法的一般做法是假设第 步的载荷 和相应的位移 为已知，然后将载荷增加为 ，如果每一步的载荷增量 足够小，则解的收敛性可以保证。
使用增量法的（1）式改写为如下形式（位移均以单自由度为例）
,其中 是用来表示载荷变化的参数，将上式对 求导，可以得到


其中 为定义的切线刚度矩阵。
（4）式是一典型的常微分方程组问题，下面介绍的是有限元中对这类方程组的求解方法
1.欧拉方法(单自由度为例)

显然为了满足精度要求， 必须是足够小的量。使用Runge-Kutte方法可以改善精度
此方法会导致解的漂移（与真实曲线上的解产生较大的误差）
为了克服每一增量步解的误差可能导致的解的漂移，可以将（5）式改写为

这里解释一下
此方法称为考虑平衡矫正的欧拉增量方法，即将前增量步的载荷和内力的不平衡量合并到当前增量步求解，一定程度上避免了解的漂移。

2.N-R方法
为了改进欧拉法的精度，现在更多采用N-R方法，如果采用N-R方法，是在每一增量步内进行迭代
则对于 的 次增量步的第 次迭代可以表示为

表示成前述的N-R形式为
采用mN-R迭代

利用mN-R方法求解非线性方程组时，可以避免每次迭代重新形成和求逆切线矩阵，但是降低了收敛速度，特别是 曲线突然趋于平坦的情况
有Aitken加速的迭代和无Aitken加速的迭代如图所示

在研究载荷增量步长自动选择的方法时，首先是假设载荷的分布模式是给定的，变化的只是他的幅值，在此情况下，外载荷可以表示为：

等效结点载荷向量可以表示为：

是载荷幅值，载荷分布实际是 的分布

求解上述载荷分布方程依然是一个N-R迭代过程，设

是经过n+1次迭代修正后得到的 数值
下面是几种常用的自动选择载荷步长的方法

此方法对于计算结构的极限载荷很有效，利用本步刚度参数可以使步长调整的比较合理，并可以减少总的增量步数，适合于计算由理想弹塑性材料组成结构的极限载荷情况。

第i增量步结构的总体刚度可用下式度量

初始结构总体刚度的度量是

其中 和 是载荷向量和按弹性分析得到的位移向量。

(1) 利用”本步刚度参数“的规定变化量自动选择增量步长。

称为第i步的”本步刚度参数“，它代表结构本身的刚度性质，与载荷增量大小无关，当结构处于完全弹性时， ，随着塑性区扩大，结构变软， 逐渐减小，到达极限载荷时，
对于比例加载的情况：

是极限载荷参数， 是 时的结点载荷向量
利用本步刚度参数可以使步长调整得比较合理，并可以减少总的增量步
(2) 增量步长的自动选择
利用(1)得到的”本步刚度参数“的规定变化量自动选择增量步长，具体的算法步骤如下
1、求弹性极限载荷参数
先施加任意载荷 ，假定结构为完全弹性求解，求出结构的最大等效应力 ，则有

是材料的初始屈服应力
2、给定第一步载荷参数增量 的值可以事先给定，用 求解第一增量步
3、给定第二及以后各增量步的刚度参数变化的预测值 ，其大小决定步长的大小，例如可在0.05~0.2之间选择，并给定刚度的最小允许值 ，则每步载荷参数增量为

然后用 求解第 增量步
然后可以通过(1)中的公式计算出”本步刚度参数“

以 迭代为例，它可以表示为

在载荷增量步长自动控制的求解方法中， 可以表示成

其中

以上两式中:

因此（7）可以改写为
其中
是在第 次迭代中由某个规定的约束条件来确定的载荷因子增量 的第 次修正量。在现在的方法中，这个条件就是某个结点的位移增量的大小，例如规定 中的最大分量 是给定的，此条件可以表示为

其中 是 的规定值， 是除第g个元素为1，其余元素为零的向量，具体迭代算法步骤如下：
(1)计算对于节点载荷模式 的位移模式 ，即

(2)计算对于不平衡结点力向量 的位移增量修正值 和 次迭代后位移增量修正值的全量 ，即


其中 是待定载荷因子增量的修正值
(3)利用条件(8)确定 ，由于

从上式可以得到

这样就确定了第 次迭代后的载荷因子增量
载荷因子增量求出之后，可以知道每一步的外载荷，从而使用 法迭代得到位移，进而求出应变，应力和当前增量步的内力等物理量。并检验收敛准则是否满足，如未满足，回到步骤二进行新的一次迭代，直至收敛准则满足为止。 关于每一个增量步某个指定结点位移增量 本身的选择，通常的方法是第1个增量步可以由某个给定的载荷因子增量 (例如令 ，其中p_e是弹性极限载荷因子， 可取5~10)，通过求解得到 ，以后各增量步的 可由下式确定，即

其中

以上是王勖成的有限单元法给出的非线性方程组解法以及一些提高运算效率的策略，如有补充或者理解偏差，请联系指正。

7. 如何解矩阵非线性方程组

如果用解线性方程组的方法求矩阵的逆，可以这样做
分别求出Ax=λi的解（其中λi表示第i个分量为1，其余分量为0的单位列向量），得到解向量xi
然后把解向量x1,x2,,xn拼接，得到的n阶矩阵就是逆矩阵。