工程師和科學家高階數學/線法

線法

簡介

線法是一種用於求解偏微分方程的有趣數值方法。其思想是將PDE半離散化為一個龐大的（連續且相互依賴的）ODE系統，然後可以使用您所熟悉和喜愛的方法，如向前推進的龍格-庫塔方法，來求解該系統。這就是“線法”名稱的由來：該解由一系列的線（更準確地說，曲線）構成。

線法僅適用於IBVP。具有邊界約束的變數（或變數）被離散化，而與初始值相關的（一個）變數成為ODE解的自變數。更正式地說，當除一個變數之外的所有變數都被離散化時，就會得到線法，從而產生一個耦合的ODE系統。

純粹的BVP，例如在某個可怕邊界上的拉普拉斯方程，可以用類似於雅可比迭代的方法求解，稱為假瞬態法。

示例問題

考慮以下非線性IBVP，其中 $h=h(z,t)$

${\frac {\partial h}{\partial t}}=2\left({\frac {\partial h}{\partial z}}\right)^{2}+h{\frac {\partial ^{2}h}{\partial z^{2}}}\,$

$h(z,0)=0\,$

$h(0,t)=1\qquad {\mbox{and}}\qquad h(1,t)=0\,$

這個令人望而生畏的 無量綱 系統描述了液體在內部角落（“凹槽”）中的流動，其中 $h$ 是液體的的高度，而 $z$ 是沿角落軸線的距離。IBVP 指出該角落最初是空的，並且在一端 ($z=0$) 的高度固定為 1。當然，液體將沿 $z$ 流動，增加 $h$ 直到到達另一邊界 ($z=1$)，其中 $h$ 始終為零。值得注意的是，雖然這個IBVP無法解析地求解，但非線性PDE的一些其他情況存在相似解，例如z=0處的固定高度，以及該位置右側無約束。

邊界值與空間變數相關聯，$z$. 假設 $h$ 在空間上離散化，但在時間上不離散化，在均勻網格上離散化為 $n$ 個點。 那麼 $h(z,t)$ 可以近似如下

$h(i\Delta z,t)\approx h_{i}(t)\,$

所以 $h_{i}(t)$ 是一個序列（或向量），它有一個索引 $i$，而不是對 $z$ 的連續依賴；但要注意，整個序列仍然對時間有連續依賴。 $i$ 是一個從 0 到 $n-1$ 的索引，注意使用了從零開始的索引，因此 $i=0$ 對應於 $z=0$，而 $i=n-1$ 對應於 $z=1$. 觀察邊界，我們知道

$h_{0}(t)=1\qquad {\mbox{and}}\qquad h_{n-1}(t)=0\,$

邊界內的序列點呢？ 最初（在 $t=0$ 時），

$h_{i}(0)=0\,$

現在假設 PDE 的右側按照您喜歡的任何方式離散化，因此

$f_{i}={\mbox{The spatial discretization of}}\quad 2\left({\frac {\partial h}{\partial z}}\right)^{2}+h{\frac {\partial ^{2}h}{\partial z^{2}}}\,}$

例如，如果使用中心差分，最終將得到

$h(z,t)$ 替換為 $h_{i}(t)$，然後用其離散化替換方程的右邊（精確的等式是對符號的輕微濫用）

${\frac {\partial }{\partial t}}(h_{i}(t))=f_{i}\,$

由於 $h_{i}(t)$ 連續地依賴於時間，而不是其他任何東西，所以這個微分方程變為

${\frac {dh_{i}}{dt}}=f_{i}\,$

將所有內容放在一起，$h_{i}(t)$ 的解是以下 IVP 的解

$h_{i}(0)=0\qquad {\mbox{(this is the initial value)}}\,$

${\begin{cases}h_{0}=1\qquad {\mbox{(this is the left boundary value)}}\\{\cfrac {dh_{1}}{dt}}=f_{1}\\\quad \vdots \\{\cfrac {dh_{i}}{dt}}=f_{i}\\\quad \vdots \\{\cfrac {dh_{n-2}}{dt}}=f_{n-2}\\h_{n-1}=0\qquad {\mbox{(this is the right boundary value)}}\end{cases}}$

解決這個問題可以得到$h(x,t)$的近似值。注意，如果使用二階前向步進法，例如二階龍格-庫塔法，該求解方法將具有與Crank-Nicolson方法大約相同的精度，但不會同時求解一堆代數方程，這將使非線性問題難以解決，因為整潔的矩陣解將無法使用。

線法在電磁學中尤為流行，例如，使用亥姆霍茲方程模擬光線穿過透鏡以獲得更好的透鏡設計；流行的原因是該系統的一階微分方程組可以（在這種情況下）進行解析求解，因此精度僅受空間離散化的限制。

需要注意的是：顯式前向步進線法解與前向步進有限差分法類似；因此，沒有理由相信線法不會出現相同的穩定性問題。

C語言中的三階TVD RK方案

以下是C語言中實現的有效TVD RK方案。記憶體會根據需要自動分配（但它只會假設記憶體分配沒有問題），並在n = 0時釋放。請注意，除非$f(x)$的離散化也是TVD（總變差遞減），否則解不會是TVD。

當然，這不是嚴格意義上的線法求解器；它可能在需要向前推進一些大型相互依賴的向量一階微分方程組時找到用途。

/*
Third order TVD autonomous RK solution. Steps x'(t) = f(x) forward, where x is some vector of size n.
 
x is the address of the first element. x must be an array of size n. f is a function as shown above.
It's fed the address of the first element of x and the index of the element being worked on, n.
Use n = 0 to free memory.
*/
void RK3(double (*f)(double *x, unsigned int i, unsigned int n), double *x, unsigned int n, double delta_t){
 	static double *x1 = NULL, *x2 = NULL;
	static unsigned int N = 0;
	if(n == 0){
		free(x1);
		free(x2);
		x1 = NULL;
		x2 = NULL;
		return;
	}
	if(n > N){
		N = n;
		x1 = realloc(x1, n*sizeof(double));
		x2 = realloc(x2, n*sizeof(double));
	}

	unsigned int i;

	for(i = 0; i != n; i++)
		x1[i] = x[i] + delta_t*f(x, i, n);

	for(i = 0; i != n; i++)
		x2[i] = 3.0/4.0*x[i] + 1.0/4.0*(x1[i] + delta_t*f(x1, i, n));

	for(i = 0; i != n; i++)
		x[i] = 1.0/3.0*x[i] + 2.0/3.0*(x2[i] + delta_t*f(x2, i, n));

}

Fortran 90等價物

subroutine RK3(f, x, n, delta_t, pass)
   double precision, external:: f
   integer         :: n
   double precision:: x(n)
   double precision:: delta_t
   double precision:: pass(*) !-- Parameters for the routine f()
   ! Locals
   double precision:: x1(n), x2(n)

   if (n .eq. 0) return
   x1 = x + delta_t *  f(x, n, pass)
   x2 = 3d0/4d0 * x + 1d0/4d0 * (x1 + delta_t * f(x1, n, pass)) 
   x  = 1d0/3d0 * x + 2d0/3d0 * (x2 + delta_t * f(x2, n, pass)) 
end subroutine RK3

f()應該定義為

function f(x, n, pass) result (r)
integer:: n
double precision:: f(n)
double precision:: pass(*)
double precision:: r(n)
! START OF CODE