方程組求解的案例
【數(shù)值算法】共軛梯度法求解線性方程組
在有限元程序開發(fā)中，線性方程組的求解是一個(gè)重要組成部分。在百萬(wàn)自由度大規(guī)模計(jì)算的情況下，線性方程組的高效快速求解對(duì)整個(gè)求解器的計(jì)算效率有著至關(guān)重要的作用。無(wú)論實(shí)際上計(jì)算的是線性問(wèn)題，還是各種非線性問(wèn)題，最終都需要落實(shí)到線性方程組的求解上去。非線性方程組的求解實(shí)際上往往就是多次求解線性方程組。
目前，線性方程組的求解主要分為直接法和迭代法兩種。
在之前的文章[數(shù)值算法與編程]高斯消去法中，我們討論的高斯消去法就是直接法的一種。而本文即將討論的共軛梯度法，是迭代法的一種，并且，其屬于目前求解對(duì)稱線性方程組的主要迭代方法。各大商業(yè)有限元軟件，在面臨對(duì)稱線性方程組的求解時(shí)幾乎都會(huì)選用各種變化形式的共軛梯度法進(jìn)行求解。
共軛梯度法的具體原理和算法如下：
假定要求解的對(duì)稱線性方程組是：
其中，A是對(duì)稱正定的系數(shù)矩陣。
則實(shí)際上待求的解也是方程
取得最小值的時(shí)候的解。
求該方程的最小值的常見方法是最速下降法，該方法算法偽代碼如下：
該方法實(shí)際上是沿著負(fù)梯度方向進(jìn)行搜索，直至殘量接近0，較為簡(jiǎn)便，但是在條件數(shù)很大時(shí)，該方法收斂很慢。
展開 大型稀疏線性方程組求解技術(shù)——工業(yè)仿真的底層核心
背景
在工業(yè)仿真領(lǐng)域，對(duì)各種現(xiàn)實(shí)世界的問(wèn)題進(jìn)行數(shù)值模擬時(shí)，如流體動(dòng)力學(xué)分析、電磁場(chǎng)仿真、結(jié)構(gòu)力學(xué)應(yīng)力應(yīng)變分析等，其控制方程通常是偏微分方程組，在經(jīng)過(guò)不同方法的隱式離散之后最終都可轉(zhuǎn)化為大型稀疏線性方程組。隨著人們對(duì)計(jì)算精度要求的不斷提高，方程組的階數(shù)也從上千階、幾十萬(wàn)階提高到百萬(wàn)、千萬(wàn)階甚至更高，所需的計(jì)算量以及存儲(chǔ)需求也隨之迅速膨脹。根據(jù)一般經(jīng)驗(yàn)，方程組求解時(shí)間會(huì)占總計(jì)算時(shí)間的70%以上，往往是整個(gè)計(jì)算過(guò)程中的性能瓶頸。如果說(shuō)求解器是工業(yè)CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術(shù)將毫無(wú)疑問(wèn)是底層求解器的核心。
NASA翼型網(wǎng)格經(jīng)過(guò)離散得到的稀疏矩陣（素材來(lái)源于網(wǎng)絡(luò)）
方法
眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優(yōu)劣，特點(diǎn)比較如下：
迭代法[1]：
對(duì)于不同類型稀疏矩陣表現(xiàn)差異較大，存在收斂性與收斂速度問(wèn)題，催生了許多預(yù)處理技術(shù)（Preconditioners）；
對(duì)原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運(yùn)算；
內(nèi)存需求小，求解速度較快，算法復(fù)雜度低；
較易實(shí)現(xiàn)并行化。
直接法[2]：
通用、穩(wěn)定；通過(guò)前后處理，能夠保證計(jì)算的收斂性與精度；
對(duì)原矩陣的編輯多（分解、排序、縮放等）；
內(nèi)存需求大，求解速度慢，算法復(fù)雜度更高；
并行度有限。
其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。
展開 技術(shù)分享︱大型稀疏線性方程組求解技術(shù)——工業(yè)仿真的底層核心
一、背景
在工業(yè)仿真領(lǐng)域，對(duì)各種現(xiàn)實(shí)世界的問(wèn)題進(jìn)行數(shù)值模擬時(shí)，如流體動(dòng)力學(xué)分析、電磁場(chǎng)仿真、結(jié)構(gòu)力學(xué)應(yīng)力應(yīng)變分析等，其控制方程通常是偏微分方程組，在經(jīng)過(guò)不同方法的隱式離散之后最終都可轉(zhuǎn)化為大型稀疏線性方程組。隨著人們對(duì)計(jì)算精度要求的不斷提高，方程組的階數(shù)也從上千階、幾十萬(wàn)階提高到百萬(wàn)、千萬(wàn)階甚至更高，所需的計(jì)算量以及存儲(chǔ)需求也隨之迅速膨脹。根據(jù)一般經(jīng)驗(yàn)，方程組求解時(shí)間會(huì)占總計(jì)算時(shí)間的70%以上，往往是整個(gè)計(jì)算過(guò)程中的性能瓶頸。如果說(shuō)求解器是工業(yè)CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術(shù)將毫無(wú)疑問(wèn)是底層求解器的核心。
NASA翼型網(wǎng)格經(jīng)過(guò)離散得到的稀疏矩陣（素材來(lái)源于網(wǎng)絡(luò)）
二、方法
眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優(yōu)劣，特點(diǎn)比較如下：
迭代法[1]：
對(duì)于不同類型稀疏矩陣表現(xiàn)差異較大，存在收斂性與收斂速度問(wèn)題，催生了許多預(yù)處理技術(shù)（Preconditioners）；
對(duì)原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運(yùn)算；
內(nèi)存需求小，求解速度較快，算法復(fù)雜度低；
較易實(shí)現(xiàn)并行化。
直接法[2]：
通用、穩(wěn)定；通過(guò)前后處理，能夠保證計(jì)算的收斂性與精度；
對(duì)原矩陣的編輯多（分解、排序、縮放等）；
內(nèi)存需求大，求解速度慢，算法復(fù)雜度更高；
并行度有限。
其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。
展開 C語(yǔ)言實(shí)現(xiàn)方程組求解 - 列主元高斯消去法和LU分解
首先是頭文件數(shù)據(jù)信息
然后是要進(jìn)行基本操作的函數(shù)文件：
==> 這里調(diào)用的上三角形方程組求解的實(shí)現(xiàn)過(guò)程如下：
主函數(shù)的測(cè)試案例為：
==> 下面進(jìn)行矩陣的LU分解，進(jìn)而求解方程組.
==> 關(guān)于下三角形方程組求解函數(shù)試下如下：
==> 這里的兩個(gè)工具類函數(shù)，如下所示：
==>
高斯消去的計(jì)算結(jié)果為：
LU分解的計(jì)算結(jié)果為：
==> 標(biāo)準(zhǔn)數(shù)據(jù)為：
標(biāo)準(zhǔn)結(jié)果為：
==> 我們所完成的兩種方法的計(jì)算結(jié)果是正確的。
支付寶贊助
微信：
【數(shù)值算法】系數(shù)矩陣非對(duì)稱時(shí)，線性方程組如何求解？-穩(wěn)定雙共軛梯度法(Bicgstab)求解線性方程組
enddo
write(*,*)"the solution of equation:"
write(*,"(es18.8)")x
end subroutine bicgstab
依據(jù)上述過(guò)程編寫程序，計(jì)算前述非對(duì)稱矩陣線性方程組求解結(jié)果：
采用matlab求解該方程組的解：
通過(guò)對(duì)比可知11次迭代已經(jīng)獲得即為準(zhǔn)確的結(jié)果。實(shí)際上，對(duì)于該方法也可以通過(guò)一定的預(yù)處理方式，使得其所需要的迭代次數(shù)更少。
以上，就是穩(wěn)定雙共軛梯度法求解線性方程組的內(nèi)容，感謝您的閱讀！
歡迎關(guān)注公眾號(hào) 有限元術(shù)
展開 大規(guī)模稀疏矩陣線性方程組求解可以有多快！
對(duì)于工業(yè)軟件研發(fā)尤其是CAE軟件研發(fā)來(lái)說(shuō)，線性方程組的求解精度和速度較為重要，在方程組規(guī)模上來(lái)以后，以Krylov子空間為基礎(chǔ)的PCG，Bicgstab，Gmres等方法相對(duì)于直接法在求解效率和精度控制上有較大優(yōu)勢(shì)。這也是abaqus等商業(yè)軟件在對(duì)于規(guī)模較大的模型時(shí)采用迭代法求解通常會(huì)比直接法求解更快的其中一個(gè)原因。
以自由度數(shù)為41w的如下分布的稀疏矩陣為例，
該稀疏矩陣形成的系數(shù)方程組在matlab中采用直接求解為1.5s,采用一定預(yù)處理下的PCG求解時(shí)間為0.55s。
用筆者自行開發(fā)的稀疏矩陣PCG求解器運(yùn)行此算例，方程組求解時(shí)間為0.628s。
此時(shí)，整個(gè)程序的運(yùn)行時(shí)間瓶頸實(shí)際上并不在方程求解，而在于從文件中讀取稀疏矩陣對(duì)應(yīng)的數(shù)據(jù)。一般情況下采用fscanf讀取數(shù)據(jù)會(huì)快于用fstream讀取。改用fscanf讀取數(shù)據(jù)后，程序總運(yùn)行時(shí)間從原來(lái)的16s變?yōu)?s。
【完】
歡迎關(guān)注公眾號(hào) 有限元術(shù)
一個(gè)講有限元技術(shù)的公眾號(hào)
展開 單元細(xì)分源代碼和線性方程組求解器
最近寫程序，需要一個(gè)高效簡(jiǎn)單的線性方程組求解器，就在別人的基礎(chǔ)上改造了一套代碼。
這套有兩個(gè)版本:fortran/C++，主要用于求解稀疏對(duì)稱矩陣；
C++ version 來(lái)源于網(wǎng)絡(luò)。 fortran version 由Esimulate改寫。 fortran version 的輸入文件為C++ version\LDL\Matrix 中任意文件，名字固定為A01。
ldlsolver.rar
C語(yǔ)言實(shí)現(xiàn)方程組求解-指針操作
在有限元的計(jì)算過(guò)程中，當(dāng)我們得出剛度矩陣之后，接下來(lái)是帶入邊界條件來(lái)求解方程組中的未知量。
1. 求解方程組中的未知量的具體方式分直接法和迭代法，一般在有限元計(jì)算過(guò)程中都采用的是迭代法，但是對(duì)于那些很小的問(wèn)題便可以采用直接法來(lái)直接得到。
2. 本次主要討論的是列主元高斯消去法。
==》 順道吐槽一下這個(gè)自帶的插入代碼的功能，很差勁。
轉(zhuǎn)貼：變參數(shù)非線性方程組的求解！
變參數(shù)非線性方程組的求解！
對(duì)于求解非線性方程組一般用fsolve命令就可以了，但是對(duì)于方程組中某一系數(shù)是變化的，該怎么求呢？
%定義方程組如下，其中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))];
%求解過(guò)程
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)
展開 如何采用simulink求解常微分方程組
通常來(lái)說(shuō)，求解一個(gè)系統(tǒng)的話采用常微分方程組去做。前面也有采用scipy進(jìn)行了常微分方程組的求解簡(jiǎn)單介紹，當(dāng)然需要用到Python。其實(shí)完全可以不用任何代碼，只用一些simulink模塊以搭積木的形式完成這個(gè)過(guò)程，而且還會(huì)方便很多。下面就介紹一下相關(guān)的方法。
所用到的核心模塊其實(shí)就是integrate模塊，只需要啟動(dòng)matlab打開simulink然后脫出一個(gè)該模塊就可以了。
首先以如下方程為例，假設(shè)初始值為0，求解區(qū)間為【0-10】
采用如下的方式搭建
simulink中的模塊求解的結(jié)果
當(dāng)然這個(gè)有點(diǎn)簡(jiǎn)單，來(lái)一個(gè)稍微復(fù)雜一點(diǎn)的
計(jì)算過(guò)程的模塊搭建如下
simulink中的模塊
計(jì)算結(jié)果如下
simulink中求解結(jié)果
當(dāng)然完全完全可以求解更加復(fù)雜的問(wèn)題，比如以下面的一個(gè)方程組為例
那么他的搭建模塊如下所示
方程組越大，則模塊會(huì)越復(fù)雜，一般可以把一部分單獨(dú)拿出來(lái)做一些封裝，然后把這個(gè)作為自己的模塊老使用，作為演示，我這里也有一個(gè)例子，就是pemfc燃料電池的例子，方程組的關(guān)系如下。
pemfc的系統(tǒng)所用到的方程
那么對(duì)應(yīng)的模塊搭建如下，可見對(duì)于較大的模型搭建還是比較難得
展開 Python 實(shí)現(xiàn)方程組的直接求解，用于求解力/位移方程組 ￥3.33
列主元三角分解法求解
千萬(wàn)自由度大規(guī)模線性方程組并行求解庫(kù)PETSc的介紹
同時(shí)，當(dāng)ksptype為preonly時(shí)，表明方程組不采用迭代求解，只采用預(yù)處理方法直接求解方程組，因此此時(shí)實(shí)際上就實(shí)現(xiàn)了線性方程組的直接法求解。
scipy求解常微分方程組
Scipy求解常微分方程組有scipy.integrate.solve_ivp和scipy.integrate.odeint，后者是較老的版本主要是采用 FORTRAN 的odepack庫(kù)里面的lsoda 方法，而前者是后面更新的函數(shù)，支持的方法也更多，按照官方的文檔介紹大致有如下的方法。
CFD學(xué)習(xí)：用時(shí)域有限差分法求解麥克斯韋方程組
要點(diǎn)
FDTD技術(shù)直接離散化麥克斯韋方程的時(shí)域偏微分形式。
頻域有限差分（FDFD）源自FDTD。
時(shí)域有限差分法是求解麥克斯韋方程組的最先進(jìn)方法，尤其是在復(fù)雜幾何形狀中。
FDTD方法可以解決與天線相關(guān)的問(wèn)題
我們經(jīng)常使用基于電流、電荷和場(chǎng)變化產(chǎn)生的電場(chǎng)和磁場(chǎng)的電器或設(shè)備。為了以數(shù)學(xué)方式表達(dá)所產(chǎn)生的電場(chǎng)和磁場(chǎng)，使用了麥克斯韋方程，并對(duì)電磁系統(tǒng)進(jìn)行了數(shù)值建模。
為了求解描述電磁場(chǎng)的方程，使用了各種數(shù)值技術(shù)。時(shí)域有限差分（FDTD）方法是解決電磁問(wèn)題最流行的技術(shù)。FDTD 方法解決了與電介質(zhì)、天線、微帶電路以及暴露于輻射的人體電磁吸收相關(guān)的問(wèn)題。在本文中，我們將深入探討 FDTD 方法。
時(shí)域有限差分 (FDTD) 方法背后的理論
FDTD方法是一種全波數(shù)值方法。該技術(shù)直接離散化麥克斯韋方程組的時(shí)域偏微分形式。為了解決電磁問(wèn)題，我們的想法是在時(shí)間和空間上使用中心差分近似來(lái)離散麥克斯韋方程組。
FDTD 技術(shù)首先由 KS Yee 通過(guò) Yee 離散方案引入計(jì)算電磁學(xué)中。在 Yee 開發(fā)的方案中，電場(chǎng)和磁場(chǎng)分量在 3 維 (3D) 空間和時(shí)間中交錯(cuò)。在所形成的3D空間中，物理電磁波傳播由法拉第定律和安培定律等值線的互連陣列來(lái)表示。使用 FDTD 技術(shù)解決電磁問(wèn)題不需要大量先驗(yàn)知識(shí)，因?yàn)?Yee 方案方法易于使用且用途廣泛。
電磁分析和 FDTD 方法
FDTD 的簡(jiǎn)單性、多功能性和靈活性使其在計(jì)算電磁應(yīng)用中廣受歡迎。由于 FDTD 方法是基于體積的，因此對(duì)于復(fù)雜結(jié)構(gòu)和介質(zhì)的建模非常有效，尤其是與有限元方法(FEM) 和矩量方法 (MOM) 相比。
FDTD 是一種時(shí)域方法。
展開 伽遼金有限元求解微分方程 -- C語(yǔ)言實(shí)現(xiàn) ￥4.5
==》 數(shù)值解解析解的對(duì)比結(jié)果圖：
==》 求解步驟主要是：
1. 寫出微分方程的弱解積分形式。
2. 進(jìn)行分布積分法。
3. 網(wǎng)格劃分。
4. 生成系數(shù)矩陣和方程組的右端項(xiàng)。
5. 進(jìn)行方程組的求解。
6. 求解出節(jié)點(diǎn)上的U值。
