牛頓下山法c語言 有限元簡單科普之 優化的牛頓法

2021-10-18 01:48:00 字數 2570 閱讀 5746

初中的時候,我們學過一篇課文,叫做《走一步,再走一步》。作者是美國的莫頓·亨特,講述了他小時候在費城,乙個酷熱的七月天裡爬懸崖的故事。小亨特爬到頂後無法下來,後來在他父親的鼓勵下通過一步一步的梯度運算,哦不,自我引導與嘗試成功地找到了下山的路徑。

很顯然,小亨特在一開始是無法找到下山路徑的解答的,只能基於初始條件,通過一步步往山下靠近,在每一步的位置座標已知的情況下才能算出下一步的向量解。

這就是有限元軟體在處理非線性模型時候的運算原理。

在我們現實中,很多事物是線性發展的,比如今年我25歲,兩年後肯定就是27歲;勻速奔跑的馬拉松運動員現在跑了20km,只要告訴我他的速度,我就能猜出十分鐘後他在什麼位置。這類線性發展的事物十分容易描述,也容易**。

但是許多事物,是非線性的,也是難以**的。比如我們人類身高隨著年齡增長的關係,程式設計師未來的收入。但是好在這類事物在短期內我們是有辦法**它的,且時間越短結果越精確。比如有人指著乙個剛出生的嬰兒問你,這孩子15年後能長到多高,你八成就答不上來;但是他問這孩子15天後能長多高,你只要調查清楚他的營養狀況和生活環境,基本能猜個**不離十。在基於彈塑性力學的有限元計算中,我們的本構模型往往是非線性的,且往往沒有簡單的方程式可以做計算。因而我們需要把計算分為很多個步長才能算出結果。

你直接問計算機在給定邊界條件下,乙個基於mcc模型的土塊直接加10kpa的力它的變形是多少,計算機一下子是無法回答你的。因為這屬於非線性計算,就好比問嬰兒十五年後的身高你也答不出來一樣。但是只要告訴計算機5kpa的時候它的變形是x(初始狀態),問5.1kpa的時候變形是多少,這時候它就可以一下子計算出來了。而且步長越小越精確,比如你問5.01kpa的時候變形多少,計算機給的結果就會精確得多。

就是靠著這麼走一步,再走一步的計算,計算機才能最終告訴你10kpa 的時候這塊土的變形是多大。因為每個步長都要經過許多次的迭代計算,所以你手裡的有限元軟體要跑好一會也是情理之中。建議大家跑有限元的話買個好的cpu,因為我曾經試過奔騰處理器和core i5的區別,速度前者大概是後者的十分之一。

廢話講得有點多,下面切入正題。

無論計算機把運算分成多少個步長(increment),每一步在本質上就是解下面這個方程:

左邊是每一步長的整體剛度矩陣乘以全域性座標中節點的位移向量增量,位移向量的下標n代表的是節點(node);右邊是荷載向量增量。角標的g表示global,右上角的i表示單個步長(increment)。

那麼這個整體剛度矩陣

解方程是我們從小學就開始學會的技能,比如一次方程,二次方程,三角函式……但是你總結一下就會發現我們基本上只能套公式去解。要是沒有求根公式呢?比如三次曲線和x軸的交點這種,那我們只能去「湊」了。我當年就是利用casio計算器裡的f(x)功能用二分法成為班裡的計算小王子(其實是太窮買不起991)。但是在計算機裡我們需要一種演算法更加簡便,能盡快能收斂到較精確解。

牛頓法,或者叫牛頓-拉夫森法(newton-raphson method)是有限元計算的一種經典方法,牛頓法的原理十分簡單,就是有乙個函式

牛頓法示意圖

牛頓法的收斂速度可謂是相當快,往往在手算中,迭代那麼三四次基本就收斂得差不多了。牛頓法的另乙個特點是精度奇高無比,每一次迭代其精度是成倍增長的——極高的精度意味著較小的累積誤差,因為「每走一步」理論上都是會有誤差的。只要是連續函式,從靠近真實解的位置出發,牛頓法基本上攻無不破。

但是這種方法也有一點小小的問題,那就是每次迭代的時候都要求一下導數。對於乙個函式的求導,在函式式複雜的情況下,是十分消耗算力的行為。所以目前的有限元軟體傾向於用一種優化的牛頓法(modified newton-raphson methods)求解:

優化的牛頓法

沒錯,就是只要第一次設定乙個斜率,然後每一次迭代的斜率都跟第一次相同。這樣雖然犧牲了收斂速度,但是免除了每步求導的過程,提高了算力。更要緊的是,這樣出錯的概率也會變小。(有些奇怪的函式求導會出問題)

通常在有限元計算中,變形很小的條件下,只要迭代那麼三到五次,解的誤差就在2%以下了。變形較大的條件下,需要迭代的次數更多。本人計算的時候一般把迭代的次數設定在500次為上限,相對誤差2%以內,這樣具有較高的運算價效比。

所以假設給一根梁加100kpa的均勻荷載,計算機會讓你先選擇步長,比如我選擇100步 (100 increments),每步1kpa,那麼每一步運算時(下圖中從a點到b點,荷載增量略誇張),梁上一點的位移和荷載曲線是這樣的:

每一步的運算(畫得略誇張,沒那麼大的)

每一步都能收斂到較高的精確值,這樣一步步位移向量增量的計算結果就算出來了。如果算到最後一步無論迭代多少次都無法收斂,原因當然是多方面的。其中乙個常見的原因就是函式在那時到了水平段了(也就是土材料fail了)。

許多傻瓜軟體都是直接幫你定好步長,然後你傻傻地坐在電腦前等結果就行了。但是其迭代運算的原理就差不多長這樣。同樣的,每一步求解都要滿足基於最小勢能原理的單元平衡方程。步長的設定對結果的精度有重要影響。

人生也是如此,走一步,再走一步。

有限元法 有限差分法 有限體積法

有限元法也叫有限單元法 finite element method,fem 是隨著電子計算機的發展而迅速發展起來的一種彈性力學問題的數值求解方法。五十年代初,它首先應用於連續體力學領域 飛機結構靜 動態特性分析中,用以求得結構的變形 應力 固有頻率以及振型。由於這種方法的有效性,有限單元法的應用已從...

有限元法,有限差分法和有限體積法的區別

原文 有限差分方法 fdm 是計算機數值模擬最早採用的方法,至今仍被廣泛運用。該方法將 求解域劃分為差分網格,用有限個網格節點代替連續的求解域。有限差分法以taylor級 數展開等方法,把控制方程中的導數用網格節點上的函式值的差商代替進行離散,從而 建立以網格節點上的值為未知數的代數方程組。該方法是...

有限元法,有限差分法和有限體積法的區別

有限差分方法 fdm 是計算機數值模擬最早採用的方法,至今仍被廣泛運用。該方法將求解域劃分為差分網格,用有限個網格節點代替連續的求解域。有限差分法以taylor級數展開等方法,把控制方程中的導數用網格節點上的函式值的差商代替進行離散,從而建立以網格節點上的值為未知數的代數方程組。該方法是一種直接將微...