線性方程組的案例
【數值算法】系數矩陣非對稱時，線性方程組如何求解？-穩定雙共軛梯度法(Bicgstab)求解線性方程組
enddo
write(*,*)"the solution of equation:"
write(*,"(es18.8)")x
end subroutine bicgstab
依據上述過程編寫程序，計算前述非對稱矩陣線性方程組求解結果：
采用matlab求解該方程組的解：
通過對比可知11次迭代已經獲得即為準確的結果。實際上，對于該方法也可以通過一定的預處理方式，使得其所需要的迭代次數更少。
以上，就是穩定雙共軛梯度法求解線性方程組的內容，感謝您的閱讀！
歡迎關注公眾號 有限元術
展開 【數值算法】共軛梯度法求解線性方程組
在有限元程序開發中，線性方程組的求解是一個重要組成部分。在百萬自由度大規模計算的情況下，線性方程組的高效快速求解對整個求解器的計算效率有著至關重要的作用。無論實際上計算的是線性問題，還是各種非線性問題，最終都需要落實到線性方程組的求解上去。非線性方程組的求解實際上往往就是多次求解線性方程組。
目前，線性方程組的求解主要分為直接法和迭代法兩種。
在之前的文章[數值算法與編程]高斯消去法中，我們討論的高斯消去法就是直接法的一種。而本文即將討論的共軛梯度法，是迭代法的一種，并且，其屬于目前求解對稱線性方程組的主要迭代方法。各大商業有限元軟件，在面臨對稱線性方程組的求解時幾乎都會選用各種變化形式的共軛梯度法進行求解。
共軛梯度法的具體原理和算法如下：
假定要求解的對稱線性方程組是：
其中，A是對稱正定的系數矩陣。
則實際上待求的解也是方程
取得最小值的時候的解。
求該方程的最小值的常見方法是最速下降法，該方法算法偽代碼如下：
該方法實際上是沿著負梯度方向進行搜索，直至殘量接近0，較為簡便，但是在條件數很大時，該方法收斂很慢。
展開 大型稀疏線性方程組求解技術——工業仿真的底層核心
背景
在工業仿真領域，對各種現實世界的問題進行數值模擬時，如流體動力學分析、電磁場仿真、結構力學應力應變分析等，其控制方程通常是偏微分方程組，在經過不同方法的隱式離散之后最終都可轉化為大型稀疏線性方程組。隨著人們對計算精度要求的不斷提高，方程組的階數也從上千階、幾十萬階提高到百萬、千萬階甚至更高，所需的計算量以及存儲需求也隨之迅速膨脹。根據一般經驗，方程組求解時間會占總計算時間的70%以上，往往是整個計算過程中的性能瓶頸。如果說求解器是工業CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術將毫無疑問是底層求解器的核心。
NASA翼型網格經過離散得到的稀疏矩陣（素材來源于網絡）
方法
眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優劣，特點比較如下：
迭代法[1]：
對于不同類型稀疏矩陣表現差異較大，存在收斂性與收斂速度問題，催生了許多預處理技術（Preconditioners）；
對原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運算；
內存需求小，求解速度較快，算法復雜度低；
較易實現并行化。
直接法[2]：
通用、穩定；通過前后處理，能夠保證計算的收斂性與精度；
對原矩陣的編輯多（分解、排序、縮放等）；
內存需求大，求解速度慢，算法復雜度更高；
并行度有限。
其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。
展開 技術分享︱大型稀疏線性方程組求解技術——工業仿真的底層核心
一、背景
在工業仿真領域，對各種現實世界的問題進行數值模擬時，如流體動力學分析、電磁場仿真、結構力學應力應變分析等，其控制方程通常是偏微分方程組，在經過不同方法的隱式離散之后最終都可轉化為大型稀疏線性方程組。隨著人們對計算精度要求的不斷提高，方程組的階數也從上千階、幾十萬階提高到百萬、千萬階甚至更高，所需的計算量以及存儲需求也隨之迅速膨脹。根據一般經驗，方程組求解時間會占總計算時間的70%以上，往往是整個計算過程中的性能瓶頸。如果說求解器是工業CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術將毫無疑問是底層求解器的核心。
NASA翼型網格經過離散得到的稀疏矩陣（素材來源于網絡）
二、方法
眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優劣，特點比較如下：
迭代法[1]：
對于不同類型稀疏矩陣表現差異較大，存在收斂性與收斂速度問題，催生了許多預處理技術（Preconditioners）；
對原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運算；
內存需求小，求解速度較快，算法復雜度低；
較易實現并行化。
直接法[2]：
通用、穩定；通過前后處理，能夠保證計算的收斂性與精度；
對原矩陣的編輯多（分解、排序、縮放等）；
內存需求大，求解速度慢，算法復雜度更高；
并行度有限。
其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。
展開 
轉貼：變參數非線性方程組的求解！
變參數非線性方程組的求解！
對于求解非線性方程組一般用fsolve命令就可以了，但是對于方程組中某一系數是變化的，該怎么求呢？
%定義方程組如下，其中k為變量
function F = myfun(x,k)
H=0.32;
Pc0=0.23;W=0.18;
F=[Pc0+H*(1+1.5*(x(1)/W-1)-0.5*(x(1)/W-1)^3)-x(2);
x(1)-k*sqrt(x(2))];
%求解過程
H=0.32;
Pc0=0.23;W=0.18;
x0 = [2*W; Pc0+2*H]; % 取初值
options = optimset('Display','off');
k=0:0.01:1; % 變量取值范圍[0 1]
for i=1:1:length(k)
kk=k(i);
x = fsolve(@(x) myfun(x,kk), x0, options);%求解非線性方程組
x1(i)=x(1);
x2(i)=x(2);
end
plot(k,x1,'-b',k,x2,'-r');
xlabel('k')
legend('x1','x2')
[ 本帖最后由 studyboy 于 2006-7-30 17:38 編輯 ]
圖片附件: k-x1.x2.bmp (2006-7-5 23:07, 689.12 K)
展開 非線性代數方程組與定理機器證明.rar
非線性代數方程組與定理機器證明.rar1理論 數學
非線性代數方程組與定理機器證明.part1.rar
非線性代數方程組與定理機器證明.part2.rar
千萬自由度大規模線性方程組并行求解庫PETSc的介紹
同時，當ksptype為preonly時，表明方程組不采用迭代求解，只采用預處理方法直接求解方程組，因此此時實際上就實現了線性方程組的直接法求解。
非線性方程組解法----討論算法實現------專題
因為要用UL迭代法計算非線性大變形問題，為提高精度，用弧長法求解非線性方程組！
多謝幫忙！
單元細分源代碼和線性方程組求解器
最近寫程序，需要一個高效簡單的線性方程組求解器，就在別人的基礎上改造了一套代碼。
這套有兩個版本:fortran/C++，主要用于求解稀疏對稱矩陣；
C++ version 來源于網絡。 fortran version 由Esimulate改寫。 fortran version 的輸入文件為C++ version\LDL\Matrix 中任意文件，名字固定為A01。
ldlsolver.rar
大規模稀疏矩陣線性方程組求解可以有多快！
對于工業軟件研發尤其是CAE軟件研發來說，線性方程組的求解精度和速度較為重要，在方程組規模上來以后，以Krylov子空間為基礎的PCG，Bicgstab，Gmres等方法相對于直接法在求解效率和精度控制上有較大優勢。這也是abaqus等商業軟件在對于規模較大的模型時采用迭代法求解通常會比直接法求解更快的其中一個原因。
以自由度數為41w的如下分布的稀疏矩陣為例，
該稀疏矩陣形成的系數方程組在matlab中采用直接求解為1.5s,采用一定預處理下的PCG求解時間為0.55s。
用筆者自行開發的稀疏矩陣PCG求解器運行此算例，方程組求解時間為0.628s。
此時，整個程序的運行時間瓶頸實際上并不在方程求解，而在于從文件中讀取稀疏矩陣對應的數據。一般情況下采用fscanf讀取數據會快于用fstream讀取。改用fscanf讀取數據后，程序總運行時間從原來的16s變為9s。
【完】
歡迎關注公眾號 有限元術
一個講有限元技術的公眾號
展開 關于Abaqus軟件求解的直接法和迭代法
針對Abaqus軟件用戶在使用軟件時，提出的在靜態隱式分析步中，方程求解的默認值為“Direct”，不是應該是“Iterative”問題，在此處做詳細的說明。
直接法：全稱為直接式線性方程求解法，該方法可以用于線性和非線性的分析，在ABAQUS/Standard模塊下，完成非線性分析時常使用牛頓方法或者其他的方法，比如弧長法，在求解的每次迭代過程中都必須要求解一系列線性方程組，而直接線性求解器就是用來尋找這些線性方程組的精確解的。ABAQUS/Standard模塊下的直接線性方程求解器使用稀疏、直接、高斯消元法，并且往往表現在分析所消耗時間的大部分時間中（尤其是大型模型的計算）—計算中方程的存儲占據著磁盤空間的最大部分。
迭代法：全稱為迭代式線性方程求解法，該方法在ABAQUS/Standard模塊下，可以用于尋找線性、非線性、準靜態、地應力、孔隙流動擴散以及熱傳導等分析步的線性方程組。由于采用迭代的技術，不能保證給定線性方程組有收斂解，當迭代求解器不收斂時，模型的改進有助于提高收斂性。在某些情況下，使用直接式線性求解可能是得到解答的唯一選擇，但當求解收斂時，使用迭代式線性求解法將獲得更精確的解答，當然這也要依賴于相對容許值的大小。通常情況下相對容許值的缺省值已經足夠精確，然而對于特殊的分析適當地調整容許值將會改善仿真的整體性能，如對于薄板或薄殼結構，相比直接式線性方程求解法，迭代式線性方程求解法將會更適合進行該結構的分析與計算。
展開 
工業軟件研發中處理超大模型(6)--有限元求解器
全文約8000字，仔細閱讀約20分鐘
在偏微分方程求解的幾種數值方法中，有限元方法(FEM）最為通用，在工程領域應用也最為廣泛，沒有之一。
FEM有時候作Finite Element Method，
有時候Finite Element Model，
注意區分。
本文將重點討論有限元方法中的隱式方法，也就是需要求解線性方程組，且線性方程組的剛度矩陣一般是稀疏對稱矩陣。還是老規矩，盡可能少用專業術語和方程，軟硬件偏軟件介紹。
將從以下幾方面討論：
1. 有限元基本原理和特點
2. 數據的生產，存儲和管理
3. 業務計算流程
4. 并行處理和分布式處理
5. 線性方程組求解方法
6. 后處理
7. 優化改進
從顯式，隱式和準靜態說起
在力學分析中常用的求解方式有兩種：顯式(Explicit)和隱式(Implicit)。在動力學分析中，尤其是對時間比較敏感的場合，普遍采用顯式方法，其特點是不需要組裝整體剛度矩陣求解線性方程組，在時間上分別對單個單元進行連續分析，時間步驟可以取得非常小，不存在收斂性問題，容易編寫程序進行并行分布式計算。
而隱式則相反，適合處理處理大跨度時間，需要組裝整體剛度矩陣，求解線性方程組，同時存在收斂性問題，但是在時間選擇上需要考慮兩步之間的特定狀況，太大誤差偏大，太小浪費計算資源。比如在熱瞬態分析中，針對兩步溫度推導方程，歐拉參數取值不同，分別使用Crank-Nicolson和Backward Euler方法。且每一步計算過后都需要動態評估時間步長是否合理。
以上圖片來自網絡
準靜態(Quasi-Static)是一種簡化的瞬態分析方法，在某些特定場合可以用較少的資源獲得不錯的計算結果。
展開 一篇文章入門“求解器”開發(全篇)
線性方程組求解
對于很多求解器來說，常用的
FEM,FVM,MOM,BEM等需要求解線性方程組。求解線性方程組一般不會自己寫，而是調用一些已有的求解庫。關于線性方程組的求解也一直是公眾號介紹的重點，參考
兩本線性方程組相關書籍pdf下載
一篇文章入門大規模線性方程組求解
大規模線性方程組解法簡介
除了以上專門介紹的線性方程組解法，在求解器相關內容介紹中，也通常會順帶介紹某些特殊的解法。
關于線性方程組求解庫的選擇，之前有過推薦：
入門Eigen
根據實際情況，可用MKL, OpenBlas，MUMPS
其它一些比如PETSc，Hypre，Trillion等稍重量級的也可以試用
某些專門加速的商業庫也可以使用
如果是需要長期從事線性方程組求解工作的，建議可以從最基礎的求解庫
BLAS,LAPACK,ScaLPAPCK用起，很多使用第三方庫碰到的問題和疑惑，根源都能從這些基礎庫里找到答案。另外很多比較老的庫也可以使用，比如SuiteSparse，SuperLu，Pardiso，TAUCS等，了解每種庫的使用特點以及適用范圍。
方程組求解的準確以及性能，和矩陣的數據特點，矩陣規模，硬件環境緊密相關。需要注意的是，在數據結構上做好封裝，方便調用不同庫之間進行切換。
在這一塊如果有繼續深入研究需求的話，一方面需要牢固的線性代數，矩陣運算基本功，另外對硬件知識，C++高性能計算也要比較熟練。高性能計算普遍在Linux環境上進行，熟悉類Linux環境是基本功。
8.
展開 大型對稱變帶寬方程組的解法
程序名稱
subroutine bandv(n,m,np,a,b,id,ir)――大型對稱變帶寬方程組的解法
功 能
本程序適用于求解大型線性代數方程組
的解，其中A為n×n階大型對稱變帶寬帶型矩陣，X和B均為n×m階矩陣，為了節省存儲量，對A采取變帶寬存儲，即每行只從第一個非零元素到對角線為止。
使用說明
子程序語句 subroutine bandv(n,m,np,a,b,id,ir)
參數說明
◎ n 整變量，輸入參數，方程組的階數。
◎ m 整變量，輸入參數，方程組右端常向量的個數。
◎ np 整變量，輸入參數，為系數矩陣A壓縮存儲的元素存儲的元素的個數。
◎ a 輸入參數，np個元素的一維數組，存放系數矩陣的壓縮存儲元素。
◎ b 輸入、輸出參數，n×m個元素的二維實數組。開始時存放方程組右端的m個n維常向量；工作結束時，存
放m個解向量，元素按列存放。
◎ id 輸入參數，n個元素的一維數組。存放系數矩陣的各個對角線元素在壓縮存儲的一維數組中的位置。
◎ ir 整變量，輸出參數，標志。若輸出值不為零，表示矩陣奇異，不求方程組的解；正常求解結束時值為
零。
方法簡介
其基本方法就是三角分解法：
只是其中的元素Kij在一維等帶寬存儲A中的位置為A（n）
分解得到L和D后，回代的求解步驟為：
◎ 由LV＝P向前回代求解V；
◎ 求對角陣D的逆陣D－1；
◎ 由V'＝ D－1V求修正分解后的載荷陣V'；
◎ 由LTX＝V'向后回代求解基本未知量X。
展開 在求解多物理場模型時，你應該選擇全耦合還是分步求解？ 附多物理場耦合模型及數值模擬導論下載
直接與迭代線性方程組求解器
無論采用全耦合方法還是分離方法，每次迭代中都會求解線性方程組。軟件為求解線性方程組提供兩類算法：
直接和迭代求解器。
直接求解器具有最穩健、最通用的優點。其缺點是需要相對大量的內存和時間，并且隨著問題規模的增加，內存需求和求解時間都會迅速增加。迭代求解器需要的內存和時間都更少，并且隨著模型大小的增加，內存需求和時間的增加速度比較緩慢。但是，迭代求解器的魯棒性較差，對于所謂的病態問題，其收斂速度較慢。例如，當材料屬性的反差非常明顯或者幾何寬高比非常大時，就會出現病態問題。近似病態的問題示例包括一根非常細長的梁的結構彎曲，或者材料電導率相差幾個數量級的電流模型。
COMSOL Multiphysics 提供的“直接”求解器包括 PARDISO、MUMPS 和 SPOOLES，以及“密集矩陣求解器”。PARDISO 或 MUMPS 的求解速度可能最快，而 SPOOLES 使用的內存可能最少。它們都應收斂到同一個解。“密集矩陣求解器”僅適用于“邊界元法”模型。
軟件中提供許多不同類型的“迭代”求解器，每個求解器都包含多個較低級別的設置。通常建議您不要手動選擇迭代求解器并調整這些設置。當已知某個特定問題適用的迭代求解器時，軟件會自動將這一組合作為選項提供。
選擇直接或迭代求解器
要在“直接”或“迭代”線性方程組求解器之間進行切換，可以轉到
全耦合特征（如果使用“全耦合”方法）或其中一個分離步驟特征（如果使用“分離”方法），并在常規欄中，將線性求解器改為其中一個可用選項。
“全耦合”特征中使用的迭代求解器。
“分離步驟”特征中使用的直接求解器。
下載地址：多物理場耦合模型及數值模擬導論
展開 