Newton的案例
基于Abaqus的Newton-Raphson算法
修正方法：在自己非線性有限元編程時，可以將上述迭代過程稍加修正，將每次迭代時的切線剛度矩陣換做初始切線剛度矩陣，并且在迭代時保持不變，即可大大減少計算量，修正方法如下圖所示：
圖 3：Newton-Raphson（N-R）迭代法與修正的 Newton-Raphson（N-R）迭代法（源自《有限元基礎教程》曾攀.P284）
求解非線性方程的Newton-Raphson算法
數學解釋
非線性方程的根記為
，按照牛頓迭代法：
牛頓迭代計算流程：
圖 4：牛頓迭代法計算流程
Newton-Raphson迭代算法：
function varargout=newton(fun,x0,ep,maxiter)
% NEWTON 牛頓法求方程的根
if nargin<4
maxiter=500; % 默認最大迭代次數
end
if nargin<3
ep=1e-8; % 默認允許誤差
end
if ~isscalar(fun)
dfun=fun{2}; % 導函數匿名函數形式
fun=fun{1}; % 函數的匿名函數形式
else
if isa(fun,'sym') % 函數以符號表達式的形式給出
dfun=matlabFunction(diff(fun)); % 導函數匿名函數形式
fun=matlabFunction
展開 非線性有限元編程 | Modified Newton–Raphson
本次分享的是基于Modified Newton–Raphson的非線性有限元編程小案例，將會介紹牛頓法的優缺點，以及如何一步步修正，最后結合彈簧單元的例子進行Modified Newton–Raphson數值實現，望大家喜歡。
Newton–Raphson有限元編程 | 一維非線性彈簧元
一維非線性彈簧模型
Newton–Raphson求解
每個彈簧元的平衡方程如下式：
整體剛度組裝，去除0左端節點項（固支條件，劃行劃列）后的整體剛度方程可表示如下：
進一步講上式化簡為1個非線性方程組
：
公式推導進行到這一步后，Newton–Raphson就可以派上用場了。有限元和數值分析緊密相關的原因就是，有限元思想將系統離散化，通過一系列的理論將物理條件轉化為一個個的數學線性/非線性方程，最終依靠數值分析的手段進行問題求解。
有關Newton–Raphson的算法細節，可參考數值分析教材或者木木前幾期的推文，本質就是一個切線逼近的過程，本次推文的目的就是帶著大家將該算法應用到具體有限元非線性求解案例中。
展開 Fortran基礎編程——牛頓拉弗森迭代
理論及在Abaqus中的實現過程已在上幾期推文基于Abaqus的Newton-Raphson算法中說明，本次主要說明Fortran編程時需要注意的地方，本文代碼主要參考：《Fortran程序設計權威指南》。
例：
主程序
主程序：Solve.f90，子程序：New_Raphson.f90，函數文件：function.f90，Module模塊：NEWTON.F90。將整個程序分塊編寫，在主程序中調用即可，大型的Fortran中顯得尤為重要。
主程序代碼：
PROGRAM main
! 主程序：Newton迭代法計算方程的根
use NEWTON
OPEN(UNIT=11,FILE='FOUT1510.TXT')
OPEN(UNIT=12,FILE='IM_RESULT1510.TXT')
CALL SOLVE(X,ITER)
WRITE(11,46)X,ITER
46 FORMAT(T5,'Newton迭代法計算方程的根',//,&
3X,'X= ',F15.10,/,&
3X,'ITER=',I5)
END PROGRAM main
由上述程序可看出，主程序只有use、call、輸出語句，甚至可以將輸出語句編一個子程序，在主程序中用call調用即可，如此以來，使得自己的代碼更加簡潔明朗，在編寫Python、Matlab也是同樣的道理。
展開 
光刻技術第20期 | 非線性壓縮感知光源-掩模優化技術及對比分析
豎直線條圖形：
（不同SMO技術對豎直線條圖形的仿真結果）
（不同SMO技術對豎直線條圖形仿真的收斂曲線）
（不同SMO技術對豎直線條圖形仿真的運行時間）
結論：
? Newton-IHTs方法PAE（3440）低于SD方法（4853）和IHTs方法（3716）。
? 收斂速度：Newton-IHTs方法收斂最快。
? 運行時間：Newton-IHTs方法效率顯著優于SD和IHTs方法。
復雜圖形：
（不同SMO技術對復雜圖形的仿真結果）
(不同SMO技術對復雜圖形的仿真收斂曲線)
(不同SMO技術對復雜圖形的仿真運行時間)
結論：
? Newton-IHTs方法成像保真度與SD方法相當，運行時間僅1161s，比SD方法（3877s）提速約3.3倍。
? 收斂速度：Newton-IHTs方法收斂最快。
? 運行時間：Newton-IHTs方法效率顯著優于SD和IHTs方法。
05/2D-DCT與2D-DFT技術的性能對比
2D-DCT基比2D-DFT基更稀疏地表示掩模圖形，可節省內存并加快計算，且對成像保真度影響較小。
(不同稀疏基上掩模圖形的稀疏系數)
(使用Newton-IHTs算法對水平條塊圖形不同稀疏基的仿真結果)
結論：
? Newton-IHTs方法的PAE（3195）遠低于SD方法（4294）和IHTs方法（3218）。
展開 有限元理論基礎及Abaqus內部實現方式研究系列4：非線性問題的求解 ￥1
2.1.2 迭代法
2.1.2.1 直接迭代法
直接迭代法取迭代方程為
2.1.2.2 Newton-Raphson法
迭代方程
詳細推導可見附錄。K的值和上面的直接法是一樣的，但最大區別是第一次迭代的u作為下一次的初始點。
2.1.2.3 Modified Newton-Raphson法
如果非線性曲線比較平緩，那么在New-Raphson迭代過程中迭代步的切線剛度矩陣可以近似為上幾個或者第一個迭代步的，即。這樣可以只在第一個迭代步計算一次K，提高每個迭代步的計算速度，但由于一般情況下
的變化較大，所以會降低收斂速度。
2.1.2.4 BFGS法
其實就是Abaqus非線性Step設置中的Quasi-Newton方法，在Abaqus做非線性可以在Newton方法和Quasi-Newton二選一。這里不做更多說明。
2.1.3 增量迭代法
一般情況增量法可以保證求解過程的收斂性但收斂速度較慢，而Newton-Raphson法收斂速度較慢但收斂性沒有保證，所以混合法結合了增量法和Newton-Raphson法來求解非線性問題。混合法首先將外載荷分成若干個增量步，在每個增量步內采用Newton-Raphson法迭代求解，在增量步內求解完成后繼續求解下一個增量步，最后將所有增量步累加起來即得到結果。
2.1.4 收斂判據
常見的收斂判據有分為兩種，失衡力準則、位移準則
2.1.4.1 失衡力準則
2.1.4.2 位移準則
2.2 Abaqus的非線性問題求解
Abaqus的非線性中上面說的三種迭代都涉及，其中General Static分析步采用Newton或者Quasi-Newton方法，而Static, Riks分析步采用Riks弧長法。我們現在只聚焦到General Static的Newton方法。
展開 深度學習中7種最優化算法的可視化與理解
, Newton_y = [], []
for it in range(1000):
Newton_x.append(x), Newton_y.append(f(x))
g = df(x)
gg = ddf(x)
x = x - g/(gg+alpha)
plt.xlim(-20, 20)
plt.ylim(-2, 10)
plt.plot(Newton_x, Newton_y, c="r", linestyle="-")
plt.scatter(Newton_x[-1],Newton_y[-1],90,marker = "x",color="g")
plt.title("Newton,alpha=%f"%(alpha))
plt.savefig("Newton,alpha=%f"%(alpha) + ".png")
plt.clf()
展開 NX NASTRAN 非線性
NX NASTRAN非線性分析提供豐富的迭代和運算控制算法，如N-R法、改進Newton法、弧長法、Newton-弧長混合法、兩點積分法、Newmark β法以及自動時間調整功能，與尺寸無關的判別準則可自動調整非平衡力、位移和能量增量，可自動完成全剛度矩陣更新或Quasi-Newton更新，使用線性搜索或自動二分載荷增量，大大提高計算效率。
NX NASTRAN可以確定屈曲和后屈曲屬性，對于屈曲問題，可以同時考慮材料和幾何非線性；對于后屈曲問題，NX NASTRAN提供三種弧長法的自適應混合使用，包括Crisfield 法、Riks法及改進的Riks法,大大提高分析效率。
?
SOL601隱式求解用于進行靜態和隱式動態高級非線性分析。適合求解低速動態問題，如金屬成型、低速碰撞分析、地震響應及生物醫學等問題。
?
SOL701顯式求解用于顯式動態高級非線性問題的分析，適合高速動態問題，如爆炸響應、高速碰撞等波的傳播問題。
SOL601與SOL701的輸入非常相似，很多情況下，用戶可以從一個分析類型的結果開始RESTAR另一個分析，滿足不同目的的數據恢復和求解，重啟動可在穩定區和非穩定區的任何一點開始。
展開 《數值方法（MATLAB版）（第四版）（英文版）——國外計算機科學教材系列》
目錄：
1 Preliminaries
1.1 Review of Calculus
1.2 Binary Numbers
1.3 Error Analysis
2 Solution of Nonlinear Equations
2.1 Iteration for Solving x=g(x)
2.2 Bracketing Methods for Locating a Root
2.3 Initial Approximation and Convergence Criteria
2.4 Newton-Raphson and Secant Methods
2.5 Aitden's Process and Steffensen's and Muller's Methods
3 Solution of Linear Systems AX=B
3.1 Introduction to Vectors and Matrices
3.2 Properties of Vectors and Matrices
3.3 Upper-Triangular Linear Systems
3.4 Gaussian Elimination and Pivoting
3.5 Triangular Factorzation
3.6 Iterative Methods for Linear Systems
3.7 Iteration for Nonlinear Systems:Seidel and Newton's Methods
4 Interpolation and Polynomial Approximation
4.1 Taylor Series and Calculation of Functions
4.2 Introduction to Interpolation
4.3 Lagrange
展開 申請兌換《數值方法（MATLAB版）（第四版）（英文版）——國外計算機科學教材系列》
Approximation
4.4 Newton Polynomials
4.5 Chebyshev Polynomials
4.6 Pade Approximations
5 Curve Fitting
6 Numerical Differentiation
7 Numerical Integration
8 Numerical Optimization
9 Solution of Differential Equations
10 Solution of Partial Differential Equations
11 Eigenvalues and Eigenvectors
Appendix:Introduction to MATLAB
Answers to Selected Exercises
Index
展開 數值方法（MATLAB版）（第四版）（英文版）——國外計算機科學教材系列
Approximation
4.4 Newton Polynomials
4.5 Chebyshev Polynomials
4.6 Pade Approximations
5 Curve Fitting
6 Numerical Differentiation
7 Numerical Integration
8 Numerical Optimization
9 Solution of Differential Equations
10 Solution of Partial Differential Equations
11 Eigenvalues and Eigenvectors
Appendix:Introduction to MATLAB
Answers to Selected Exercises
Index
展開 
ANSYS Workbench包含尖角的接觸分析案例詳解
求解后，不收斂：
Newton-Raphson殘余力：
加密網格，對參與接觸的6條邊設置edge sizing，單元尺寸0.16mm，
再次求解后，456步后收斂：
變形動畫：
===================分割線，以上為正文======================
二次開發調試技術
一般調試格式
（1）求解（solution）調試格式
執行命令：
/DEBUG,-1,F1,F2,F3,F4,F5,F6,F7,F8,F9
參數說明：
F1—1 (輸出基本求解結果控制調試結果)
F2—1 (輸出使用Newmark常數的瞬態計算調試結果)
2 (輸出使用速度與加速度的瞬態計算調試結果)
F3—1 (輸出單元矩陣調試結果，包括矩陣與載荷矢量)
2 (輸出單元矩陣調試結果，只包括載荷矢量)
3 (輸出單元矩陣調試結果，包括矩陣對角元素和載荷矢量)
F4—1 (輸出自動時間步長(auto time stepping)調試結果)
F5—1 (輸出多物理場調試結果)
F6—1 (輸出弧長(arc-length)調試結果)
F7—1 (輸出基本Newton-Raphson調試結果)
2 (輸出Newton-Raphson調試結果，包括非平衡力或增量位移或每個DOF)
3 (輸出Newton-Raphson調試結果，包括施加載荷與每個DOF 上的n-r恢復力)
F8—1,2 (輸出位移矢量以及位移指針調試結果)
2 (輸出位移矢量以及增量位移調試結果)
3 (輸出位移矢量以及接觸數據庫調試結果)
F9—1 (輸出臨時程序員調試結果)
（2）單元(element)調試格式
執行命令：
/DEBUG,-3,G1,G2,G3,G4,G5,G6,G7,G8,G9
參數說明：
G1—1 (輸出基本單元通過(pass)調試結果)
G2—1 (輸出單元位移和坐標調試結果)
G3—1 (輸出單元矩陣調試結果，包括矩陣與載荷矢量)
2 (輸出單元矩陣調試結果，只包括載荷矢量)
3 (輸出單元矩陣調試結果，包括矩陣對角元素和載荷矢量)
G4—1 (輸出單元載荷調試結果)
G5—1 (輸出單元實常數(real constant)調試結果
展開 ABAQUS/Standard vs. 非線性 (基礎篇)
[Other 選項卡]
Method：Solver求解方法，直接法和迭代法
Solution Technique: 求解技術，分為Full
Newton和Quasi-Newton， 兩者之間的最大差距就是，Full Newton在每個迭代步中結構的剛度矩陣都會重新生成，而Quasi-Newton默認下，每8個增量步才會重新生成剛度矩陣。
Extrapolation： 外推。默認Linear，非線性分析步的第一個增量步之后的每一個增量上，ABAQUS/Standard通過從前面的增量外推解來評估當前增量的解。默認情況下，使用100%的線性外推。如果Ti ≤ DE * Ti-1，則放棄外推，DE默認為0.1。
【PS】如喜歡本文請掃一下二維碼，關注一下公眾號謝謝。
展開 13篇近代數值分析經典文獻
update formula
Broyden/Fletcher/Goldfarb/Shanno (`70) BFGS quasi-Newton update formula
Karmarkar (1984) interior pt methods for linear prog.
