不求導數矩陣的非線性方程組演算法

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階矩陣就是逆矩陣。