工程师和科学家高级数学/线法

线法

简介

线法是一种用于求解偏微分方程的有趣数值方法。其思想是将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