數值算法的案例
【數值算法】共軛梯度法（二）-預處理共軛梯度法
求得系數矩陣A的不完備的下三角矩陣L后，
令
再采用預處理共軛梯度法的具體算法獲得最終解。當然，由于L是下三角矩陣，因此預處理方程一般通過兩次“回代”（參考本公眾號文章[數值算法與編程]高斯消去法的回代部分）即可求解。
以之前的文章共軛梯度法中的原始方程求解為例，采用不完備Cholesky分解預處理求解如下：
僅需3次迭代，即獲得收斂解，而原來的常規共軛梯度法需要9次迭代。并且，matlab計算結果如下：
通過對比不難看出，雖然僅僅是3次迭代，但是已經具備較高的求解精度。
在實際開發中，共軛梯度法還有較多的發揮空間，比如，比如，知名有限元大師Thomas J.R. Hughes在1983年創立了一種基于共軛梯度法的element by element算法，這種方法不需要組裝整體剛度矩陣，而是通過逐個單元進行求解，求解效率很高，且由于不需要組裝整體剛度矩陣，計算過程中的內存需求顯著減少，并且免去了常規的采用稀疏矩陣存儲有限元剛度矩陣的組裝過程，實際上，相較于常規的矩陣，對于CSR，CSC等格式的有限元整體剛度稀疏矩陣組裝，并不是一件十分容易的事情。以上，就是共軛梯度法（二）之預處理共軛梯度法的全部內容，感謝您的閱讀！
【完】
==================================================================================
參考文獻：
【1】《矩陣計算》，Gene H Golub，Charles F.
展開 【JY】淺析各動力求解算法及其算法數值阻尼(人工阻尼)
【中心差分法】由于中心差分法所需要的時間步長比較短，實質上會讓該算法的譜半徑的模長等于1，也就是該算法并不能調整數值阻尼。
[數值算法與編程]選主元高斯消去法
在之前的文章[數值算法與編程]高斯消去法 中，本公眾號編寫了高斯消去法求解線性方程組的具體代碼。其具體算法如下：
（1）消元部分
（2）回代部分
很明顯，對于某些矩陣，使用上述算法可能會出現a(i,i)為0的情況，而一旦出現這種情況，該算法實際上就無法繼續進行求解。
以以下方程組為例：
上述方程組的系數矩陣為：
第一次消元后，系數矩陣變為
顯然，由于A1(2,2)=0，下一次消元已經無法進行，因此直接采用高斯消去法對該方程組是無法進行求解的。
針對該問題，改進的算法叫選主元高斯消去法。
具體步驟如下：
1.設置增廣矩陣AB=[A,B]
2.對增廣矩陣的第i（i=1~N-1）列進行以下處理：
2.1 設amax=ABS(AB(i,i))，IDmax=i
2.2對第i列的對角線以下的元素進行遍歷，如果有元素的絕對值大于對角線元素的絕對值（amax），則獲得該元素所在行并將其行數賦值給IDMAX，并將AB矩陣中該行（即IDMAX行）與第K行進行元素交換。
對變換后的矩陣按照普通的高斯消去法進行求解。
展開 【數值算法】共軛梯度法求解線性方程組
在之前的文章[數值算法與編程]高斯消去法中，我們討論的高斯消去法就是直接法的一種。而本文即將討論的共軛梯度法，是迭代法的一種，并且，其屬于目前求解對稱線性方程組的主要迭代方法。各大商業有限元軟件，在面臨對稱線性方程組的求解時幾乎都會選用各種變化形式的共軛梯度法進行求解。
共軛梯度法的具體原理和算法如下：
假定要求解的對稱線性方程組是：
其中，A是對稱正定的系數矩陣。
則實際上待求的解也是方程
取得最小值的時候的解。
求該方程的最小值的常見方法是最速下降法，該方法算法偽代碼如下：
該方法實際上是沿著負梯度方向進行搜索，直至殘量接近0，較為簡便，但是在條件數很大時，該方法收斂很慢。
展開 
近場動力學快速入門程序——桿和板，鍵型本構及兩種求解器（顯示求解和隱式求解） ￥150
二、數值算法簡介
PD理論很難求得解析解，所以求解PD基本都是用數值算法。目前求解PD的數值算法可以分為兩類，即顯示算法和隱式算法。此外，結合離散方式的不同，求解PD的數值算法可以進行如下的劃分。
（1）無網格方法或離散粒子法
該種方法是Silling博士于2005年發表的一篇文章中提出的方法。該種方法將連續的物體離散為許多規則的有體積的質量塊，每一個質量塊都將質量集中到小塊體積的幾何中心處，那么待求解的結構就被離散為粒子系統，近場域的積分項自然而然地離散為求和的形式。當完成離散后，就可選擇是用顯示方法進行求解還是隱式方法進行求解。
對于動力學問題，常用的顯示方法是中心差分法，而對于準靜態問題，則一般使用自適應動力松弛法。自適應動力松弛法是Madenci教授于2010年發表的一篇文章中正式推廣的一種依然采用中心差分格式的方法。
對于動力學問題，常用的隱式方法是Newmark法，而對于準靜態問題，隱式方法的核心思想則是牛頓迭代法，該方法最核心的地方在于如何求解切線剛度矩陣。總的來說剛度矩陣的獲得有解析法和計算法兩種。解析法可以參考Silling博士2010年發表的論文，里面提出了模量態的概念，但由于解析法普適性沒有計算法好，所以大多數支持無網格隱式求解的開源軟件都采用計算法。計算法的思想是給一個非常小的擾動位移，然后用中心差分或向前差分或向后差分代替求導數，從而得到切線剛度矩陣。
支持離散粒子法的開源軟件有Peridgm、LAMMPS等，非開源的有Silling本人開發的EMU。
（2）有限單元法
有限單元這個詞可以形容一種理論模型，但其本身也是一種求解偏微分方程的數值方法。PD理論的平衡方程不是偏微分方程而是積分方程，所以在PD理論最開始時幾乎沒人用這種數值方法求解。
展開 北京電影學院發了一篇滿是數學公式的計算機頂會論文，并開源了其代碼（轉載）
他們的論文研究的算法，甚至直接和工業光魔合作，使用在當年星球大戰特效鏡頭的制作上
還有很多科研機構也是不斷地在算法上推陳出新，甚至像迪士尼這樣家里有礦的影視工作室還親自上陣，不斷地革新著美利堅的影視制作技術，帶動起了好萊塢產、研結合的影視工具生態。
而不管是由工業光魔發表的用來補充煙火高精度細節的算法：
還是來源于迪士尼研究院-ETH Zurich的小波湍流論文：
都是不斷在嘗試攻克一個令人頭疼的問題：如何在使用高效穩定的流體特效仿真算法前提下，盡可能地保證流體運動的細節，這最終落在了對于N-S方程中的對流部分的求解上。
在圖形學特別是影視特效制作軟件中，為了追求大時間步長的穩定性，對流方程的求解多采用半拉格朗日方法或其變種。這類方法雖然高效穩定，但有它最大的問題：數值粘性。
所謂數值粘性，顧名思義，就是物理系統中本不存在的粘性，是被人類設計的數值算法在求解方程的過程中帶入系統的，與真實的物理事實相左。我們先領略一下“數值粘性” 會帶來的問題：
比如我要沿著紅線所示的速度和方向“對流”白色的方塊
傳統的對流算法由于數值粘性會使物質模糊掉
我們的算法盡大可能地避免了數值粘性， 保留了清晰地邊界
流體的運動更為復雜，計算機算法保持清晰（守恒）的能力，就顯得尤為重要。
這是一個由傳統算法計算得到的煙霧動畫，由于數值粘性的存在，煙霧的形態被過度地模糊了。
為了對抗這種過度平滑，科技工作者們提出了混合的粒子-網格方法，然而這類方法存在著計算速度慢，計算結果差強人意的缺陷：往往由于粒子的非連續性，流體的湍流運動會在粒子間拉出空隙，最終在視覺上產生噪聲：
我們提出的算法既能最大可能地保持流體的湍流細節，又能維持流體場的連續性，有效地解決特效解算中這一老大難問題。
展開 [數值算法與編程]高斯消去法
對上述過程進行算法總結：
（1）
上述過程完成的是對左側的系數矩陣和右側的向量進行變換，n是系數矩陣a的階數，c是用于行變換的系數，通過上述變換后，得到的a為上三角矩陣，b為對應的右側向量。
（2）
上述過程完成的是依次將求得的值回代入方程中，最終獲得的b，即為方程ax=b的解。
編程：
上述例子是針對的滿秩系數矩陣的求解，200階的系數矩陣求解時間約為0.01s的級別。在實際的有限元計算中，形成的有限元剛度矩陣通常為具有對稱性的帶狀稀疏矩陣（即非0元素主要集中在矩陣對角線附近），因此針對帶狀稀疏矩陣，還可進一步對高斯消去法進行優化，使得其求解效率更高。
參考文獻：《有限元法（第二版）.下.理論、格式與求解方法》.Klaus-Jurgen Bathe著 軒建平 譯
《Introduction to finite elements in engineering-Pearson (2012)》Belegundu, Ashok D._ Chandrupatla, Tirupathi R.
歡迎關注微信公眾號 有限元術 分享有限元理論與實踐知識
展開 VOF算法的浮體入水過程的數值模擬 ￥499
采用VOF 法求解氣液界面,結合k-e湍流模型和動網格技術模擬物體入水這一簡單氣液固多相流流動。
浮體首先自由下落，流體浮力和粘性阻力在接觸液體后逐漸增大，加速度隨之減小，速度增加變緩，當浮力與粘性阻力之和等于重力時，加速度等于零，達到最大下降速度，之后開始減小，直至減小到零，達到最大入水深度；接著物體緩慢上升，粘性阻力改變方向，當再次接觸到液面后，浮力減小，速度增加變緩，當浮力等于重力與粘性阻力之和時，達到最大上升速度，之后開始減小，直至減小到零，達到最大上升高度。浮體重復上述運動過程，且由于流體粘性，幅度逐漸衰減。
基于經驗公式的不同硬度下橡膠Mooney?Rivlin模型本構參數的確定方法（使用LS-DYNA隱式算法進行準靜態橡膠壓縮數值模擬） ￥12.86
基于經驗公式的不同硬度下橡膠Mooney?Rivlin模型本構參數的確定方法
—使用LS-DYNA隱式算法進行準靜態橡膠壓縮數值模擬
一、引言
橡膠材料的力學特性一般是通過材料力學性能試驗得到應力-應變數據，之后擬合相應的本構模型來得到其材料系數，然而這組系數只能在橡膠相應的實驗應變范圍內使用，一旦超出實驗應變范圍，這組系數就不再可靠。考慮到實驗的成本、實驗條件的多變、橡膠的材料不均勻及仿真研究時的迅速、高效性，本文基于理論分析和實驗經驗結果，結合仿真分析在不需進行試驗的前提下對不同硬度的橡膠Mooney?Rivlin模型本構參數予以確定，所確定的本構參數可滿足大部分仿真工況。
Mooney?Rivlin是一個比較經典的橡膠本構模型，使用它幾乎可以模擬所有橡膠材料的力學行為，其適用于中、小變形，一般可應用于應變約為100%(拉)和30%（壓）的情況。在仿真分析中使用較簡單、應用最廣泛、精度可接受的應變能密度函數首選Mooney?Rivlin模型，其是可表達接近不可壓縮天然橡膠應力應變特性的較合理的橡膠本構模型。
二、理論分析
橡膠的剪切模量和彈性模量主要取決于其邵氏硬度，根據彈性理論：
由式（1）和（2），令彈性模量相等可得：
由于橡膠的容積彈性模數K≈2720N/mm2，剪切模量G≤2.4N/mm2,代入可得其泊松比典型值為0.4996，與0.5十分接近，本構模型參數確定時可將泊松比視為0.5。因此橡膠材料的彈性模量和剪切模量有如下關系：
Mooney?Rivlin模型的表達式為：
該模型可很好的描述變形小于150%的橡膠材料力學性能，完全能夠滿足橡膠實際應用的性能計算。
展開 國內鋰電池CAE仿真軟件的突破口
國內鋰電池CAE仿真軟件的突破口
一方面電化學建模逐步深入到電極介觀尺度，另一方面一些相對新興的數值模擬技術被用來探究介觀尺度上的電化學反應過程機理。國內鋰電池CAE軟件應當聚焦這兩個方面進行突破。同時要使用恰當的模型降階、時間離散等計算加速方法，兼顧模型的精度和計算效率。
一是不斷優化仿真模型及其控制方程，及時跟蹤電池仿真在微觀、介觀尺度上的前沿進展，將精度更好的理論模型內置在軟件之中，并開發出適配電池領域的網格劃分工具。其次是可以向COMSOL學習，保持軟件的開放性，這一點是指在建模階段仿真工程師可以自定義修改設立控制方程組（偏微分方程組）及其假設條件，而不是只能采用軟件內置的方程組。軟件不斷跟進業界的先進實踐經驗，通過加強與業界人員的合作來優化迭代自身軟件求解器的性能。相較于其他成熟的仿真領域，這一點對于電池領域的仿真軟件更為凸顯。
求解方法方面注重運用新的求解方法，并將多種數值算法結合起來，比如LBM、DEM、FEM和FVM 結合起來，并注重數值計算方法和工程實際業務場景深度結合。同時工業仿真軟件需要權衡精度和效率，這就需要恰當地使用降階模型、伽遼金投影法、時間離散方法如龍格－庫塔法（Runge-Kutta methods）等計算方法，盡可能確保在不損失仿真精度的情況下提高計算效率。
除了理論模型和數值算法，仿真軟件還需要注重工程經驗的積累，特別是對于電池這一非線性的復雜系統來說。因此，國內鋰電池CAE軟件公司一方面需要在產品功能上加強與實驗測試數據的結合，為設計和仿真人員提供更便捷的仿真服務。這里具體來說，電池仿真通常需要大量參數輸入，部分參數需要標定和實驗以及文獻參考得來，如果積累了大量的測試數據，就可以提供一個豐富的數據庫，更加方便快捷提供模擬仿真所需參數，避免因為仿真人員標定出現問題帶來較大的偏差，進一步縮短設計驗證周期。
展開 運用S-ALE算法實現殺爆戰斗部破片拋射數值模擬 ￥150
該案例基于ls-dyna R11講解S-ALE算法及運用。
包含了：
算法介紹；
算法關鍵字設置及含義；
邊界條件施加；
殺爆戰斗部破片拋射案例；
關鍵技術要點。
【數值算法】系數矩陣非對稱時，線性方程組如何求解？-穩定雙共軛梯度法(Bicgstab)求解線性方程組
在前面的文章和中表明共軛梯度法是求解對稱正定線性方程組的一種有效方法，當針對不同的系數矩陣采用不同的預處理方式時，其可以以較少的迭代次數獲得較高精度的解。然而，該方法的一個缺點就是其只能適用于對稱正定系數矩陣，當系數矩陣不再是對稱正定時，此方法可能失效。
以下舉例：
上面矩陣A為非對稱矩陣，采用共軛梯度法求解過程如下：
該方程組采用共軛梯度法迭代4862次依然未收斂。因此，對于該非對稱方程，可以認為，共軛梯度法幾乎是失效的。
在實際工程中，有限元方法形成的剛度系數以對稱正定居多，但是實際上也存在非對稱的可能，例如，當材料本構采用摩爾-庫倫本構時，其形成的剛度矩陣就有可能會是非對稱的，此時如果是使用商業軟件，應當在軟件中選擇非對稱求解器。如果是自主編程且采用迭代法求解線性方程組，則需要找到一種適用于非對稱矩陣的求解方法。
常見的非對稱系數矩陣求解方法主要有：廣義最小殘差法(GMRES)，雙共軛梯度法（Bicg）穩定雙共軛梯度法(BiCGStab)，穩定混合雙共軛梯度法（BiCGStab(l)），這些方法相對于常規的共軛梯度法在推導上均增加了一些難度，實際推導往往較為復雜。本文不展開推導，僅對穩定雙共軛梯度法(BiCGStab)的偽代碼作簡要粘貼。
BiCGStab法的具體計算過程如下：
具體代碼：
program bicgstab_main
implicit none
integer,parameter::n=8
real(8)::a(n,n),b(n)
real(8)::x(n),x0(n)
integer::i,j
integer,parameter::m=20
real(8)::c(m,m),d(
展開 21GA-ELM，遺傳算法優化ELM預測，并和優化前后以及真實數值進行對比，確定結果 ￥18.8
GA-ELM，遺傳算法優化ELM預測，并和優化前后以及真實數值進行對比，確定結果，基于MATLAB平臺，程序已經調通，可以直接運行，需要直接拍下。
數值方法(MATLAB版第4版)/國外計算機科學教材系列
數值方法(MATLAB版第4版)/國外計算機科學教材系列
作者：（美）馬修斯
市場價：￥49.00元
卓越價：￥38.80元折扣：79折節省：10.20元
VIP 價：￥38.10元折扣：77折節省：10.90元
華北地區1-2天發貨 華東地區1-2天發貨 華南地區1-2天發貨
具體確認您屬于哪個地區請點擊
用戶評分 暫無用戶評分
編輯推薦
本書全面介紹了數值方法的理論和實踐知識，注重對利用MATLAB軟件實現各種數值算法的實際能力的培養，有助于加強學生的數學理論基礎，培養學生實際處理數值計算問題的能力。書中內容豐富、覆蓋范圍廣，對于不同學習對象和學習目的，可以選擇相應的章節，形成理論與實踐相結合的學習策略。本書包含數值方法的眾多研究領域，可滿足不同專業和不同層次的學生的需求，尤其適用于數學、計算機、物理和工程專業的人員。
以實際例題清晰而深入淺出地說明概念、解釋定理；
包含大量的習題和編程題，范圍涉及多個不同的應用領域；
強調利用MATLAB進行數值方法的程序設計，包含可直接使用的代碼實例。
本書介紹了數值方法的理論及實用知識，并講述了如何利用MATLAB軟件實現各種數值算法，以便為讀者今后的學習打下堅實的數值分析與科學計算基礎。本書內容豐富，教師可以根據不同的學習對象和學習目的選擇相應的章節，形成理論與實踐相結合的學習策略。書中的每個概念均以實例說明，同時還包含大量的習題，范圍涉及多個不同領域。通過這些實例進一步說明數值方法的實際應用。本書的突出特點是強調利用MATLAB進行數值方法的程序設計，可提高讀者的實踐能力并加深對數值方法理論的理解；同時它的覆蓋范圍廣，包含數據方法的眾多研究領域，可以滿足不同專業和不同層次學生的需求。
展開 21GA-ELM，遺傳算法優化ELM預測，并和優化前后以及真實數值進行對比，確定結果。 ￥15.9
GA-ELM，遺傳算法優化ELM預測，并和優化前后以及真實數值進行對比，確定結果。基于MATLAB平臺，程序已經調通，可以直接運行。
