當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

四阶龙格库塔法解一维扩散方程

發布時間：2023/12/31 编程问答 33 豆豆

生活随笔收集整理的這篇文章主要介紹了四阶龙格库塔法解一维扩散方程小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

龍格庫塔法是解ode的利器，最近有個project是用龍格庫塔法解pde，做的時候網上這方面能搜到的資料還不多，傳一下自己的解法。

核心思路

直線法（method of line）把pde化為ode，離散化為空間n個格點上的一元函數，用向量表示

u(x,ti)→{u1(t),u2(t)...,un(t)}iT=u?iu(x,t_i)→\{u_1(t),u_2(t)...,u_n(t)\}_i^T=\vec{u}_i

空間的二階導用中心差分格式處理，寫成矩陣形式

?u2?2x=1Δx2(?2100?01?210??01?21???????1000?1?2)?{u1(t)u2(t)??un(t)}=?u?t=C?u?i\frac{\partial u^2}{\partial^2 x}=\frac{1}{\Delta x^2} \begin{pmatrix} -2 & 1 & 0 & 0 &\cdots & 0 \\ 1 & -2 & 1 & 0 &\cdots & \vdots \\ 0 & 1 & -2 & 1 & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \cdots & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{pmatrix}·\left \{ \begin{matrix} u_1(t)\\u_2(t)\\\vdots \\\vdots \\u_n(t) \end{matrix} \right\} =\frac{\partial u}{\partial t} =C·\vec{u}_i

擴散方程的時間一階偏導用4階龍格庫塔方法轉化為線性方程
四階龍格庫塔常用格式

?u?t=?2u?x2→ui+1?uih=16(k1+k2+k3+k4)=C?ui?\frac{\partial u}{\partial t} = \frac{\partial^2u}{\partial x^2}\to \frac{u_{i+1}-u_{i}}{h}=\frac{1}{6}(k_1+k_2+k_3+k_4) =C·\vec{u_i}

4.中間量取值

這是最關鍵的一步。龍格庫塔解ode之所以方便，因為f是表達式已知的顯函數，直接代入完事。但pde里面偏導表達式如果有還好辦，沒有的話像擴散方程得換條思路。注意到k是時間偏導的函數值，取決于x和y。在擴散方程中k其實是空間二階偏導的在某一點的函數值。雖然定義上某一點y的二階偏導由x,y共同決定，但實際上把二階偏導視為微分算符的話，偏導的值由y便可確定。換言之，擴散方程中u是一個隱函數，不同的k對應不同的u。\begin{matrix} {\color{Red} 這是最關鍵的一步。}\\ {\color{Red} 龍格庫塔解ode之所以方便，因為f是表達式已知的顯函數，直接代入完事。}\\ {\color{Red}但pde里面偏導表達式如果有還好辦，沒有的話像擴散方程得換條思路。} \\ {\color{Red}注意到k是時間偏導的函數值，取決于x和y。} \\ {\color{Red}在擴散方程中k其實是空間二階偏導的在某一點的函數值。} \\ {\color{Red}雖然定義上某一點y的二階偏導由x,y共同決定，}\\ {\color{Red}但實際上把二階偏導視為微分算符的話，偏導的值由y便可確定。}\\ {\color{Red}換言之，擴散方程中u是一個隱函數，不同的k對應不同的u。}\\ \end{matrix}

$k1=C?ui?,k2=C?（ui?+hk12）,k3=C?（ui?+hk22）,k4=C?（ui?+hk32）k_1=C·\vec{u_i},k_2=C·（\vec{u_i}+\frac{hk_1}{2} ）,k_3=C·（\vec{u_i}+\frac{hk_2}{2} ）,k_4=C·（\vec{u_i}+\frac{hk_3}{2} ）$

迭代
至此，我們得到一個在時間方向上迭代的前向差分格式

ui+1=ui+h6(k1+k2+k3+k4){u_{i+1}=u_{i}}+\frac{h}{6}(k_1+k_2+k_3+k_4)

C++用一個循環實現，調用eigen線代庫進行矩陣運算，vector容器定義循環的列向量，ofstream把每次循環結果列向量的值裝進excel表。

計算結果

邊界與初始條件，以及計算范圍與步長
$D=0.1\\ u(0,t)=10^{13},u(x>0,0)=0 \\ t∈[0,5],tsteps=10000;x∈[0,1],xsteps=50$
小步長多步計算，基本等于恒定源擴散的解析解——余誤差函數

縮小范圍，減少步數，計算誤差
$D=0.01t∈[0,3],tsteps=300;x∈[0,1],xsteps=30D=0.01\\ t∈[0,3],tsteps=300;x∈[0,1],xsteps=30$
龍格庫塔算出來的誤差隨函數值增大而增大，空間越往后誤差越大，時間越往后誤差越小

對比與歐拉法的誤差

誤差龍格庫塔-歐拉<0，說明確實比歐拉法精確。

總結

以上是生活随笔為你收集整理的四阶龙格库塔法解一维扩散方程的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。