计算科学 - 这个扩散方程组是适定的吗？ - 吾爱随笔录

这个扩散方程组是适定的吗？

计算科学 pde 扩散曲柄尼科尔森体贴入微

2021-11-29 19:43:04

我正在使用标准的 Crank-Nicholson 算法来求解这个由两个耦合扩散方程组成的系统：

\begin{matrix} (1) & \dot{u} - \dot{v} = \frac{\partial}{\partial x} (α (x) \frac{\partial u}{\partial x}) \end{matrix}

$\dot{u} - \dot{v} = \frac{\partial}{\partial x} \left( \alpha(x) \frac{\partial u}{\partial x} \right) \label{1} \tag{1}$

\begin{matrix} (2) & \dot{v} - \dot{u} = \frac{\partial}{\partial x} (β (x) \frac{\partial v}{\partial x}) \end{matrix}

$\dot{v} - \dot{u} = \frac{\partial}{\partial x} \left( \beta(x) \frac{\partial v}{\partial x} \right) \label{2} \tag{2}$

在哪里 $\alpha$ 和 $\beta$ 是完全确定的函数 $x$ . BC 只是狄利克雷条件。我将解向量排序为 $u_1, v_1, u_2, v_2, \ldots$ 得到一个带矩阵。

线性系统通过写离散化得到 $\eqref{1}$ 然后 $\eqref{2}$ 对于第一个节点，然后是第二个节点的这两个方程，依此类推。

我得到的结果并非完全没有意义，但振荡如此之大，以至于无法利用数据。基本上，解决方案在一半的时间步中会相当不错，而其他的将是极低的负值，几乎没有物理意义。换的时候没有改善 $\Delta t$ 或者 $\Delta x$ 或尝试不同的 BC。

我无法理解这种现象的来源。这个问题可能是不适定的吗？还是与方程式的编写/排序方式有关？

2个回答

我认为它可能是不适定的，因为与时间相关的部分是线性相关的。如果将两个与时间相关的方程加在一起，则会得到一个与时间无关的方程：

$(\alpha(x)u_x)_x + (\beta(x) v_x)_x = 0.$

你的初始条件是否满足这个方程？

我对你的问题很好奇，对我来说测试你的两个方程的 CN 方法并不难。添加 - 以“添加”开头的部分是为了解决评论而写的。

大卫给出的观点很重要。这是初始条件的三个图（添加 - 不满足建议的条件 $\alpha u_{xx} + \beta v_{xx}=0$ )，第一个时间步和第五个时间步，用隐式欧拉法计算，蓝色曲线是u，另一个是v：

扩散系数 $\alpha=\beta=1$ 是恒定的。添加 - 对于不寻常的任意小时间步，解决方案在第一个时间步中快速变化。之后可以观察到标准数值扩散（第 5 个时间步）。

正如您所观察到的（第 1 和第 5 时间步），使用 Crank-Nicolson 进行相同的计算很奇怪：

我想有人可以稍后直接对您的 PDE 发表评论，但是从这个简化的数值分析中，我猜您没有机会使用 Crank-Nicolson 方法获得合理的解决方案（添加 - 如果您的初始功能不满足 $(\alpha u_x)_x + (\beta v_x)_x=0$ ）。我怀疑使用隐式欧拉方法的数值解是否合理，但如果是，则表明对于不满足大卫描述的此 PDE 的初始条件，您在解中获得任意小时间的跳跃，因此您有一些“无限”的进化速度。已添加 - 如注释中所述，对于显式 Euler 方法，情况甚至更糟。

这只是启发式答案，但至少它证实了您对 Crank-Nicolson 方法对您的问题的怀疑。下面是 Mathematica 中生成图片的代码。

添加 - 如果您的初始函数满足所需的固定 PDE（但使用例如提供的代码中的有限差分法的近似形式），那么您可以使用隐式 ( $\theta=1$ 在代码中）或显式（ $\theta=0$ ) Euler 方法或 Crank-Nicolson 方法 ( $\theta=0.5$ ) 对于时间离散化，它们将起作用，因为如果初始函数满足演化数值函数，这些方法也保留此属性。

Model data
t0 = 0.; T = 0.1; L = 10.; alpha = 1.; beta = 1.; 
u0[x_] := If[Abs[x - 3] <= 1., 1., 0.]; 
v0[x_] := If[Abs[x - 7] <= 1., -1., 0.]; 

Method's input
Nd = 5; M = 20; theta = 0.5; \[Tau] = T/Nd; h = L/(M + 1); 
Mat = SparseArray[{}, {2*M, 2*M}]; 
ps = Array[{}, {2*M}]; 

Method
For[i = 1, i <= M, i++, If[i != 1, Mat[[i,i - 1]] = (-alpha)*\[Tau]*theta]; Mat[[i,i]] = h^2 + 2*alpha*\[Tau]*theta; 
   If[i != M, Mat[[i,i + 1]] = (-alpha)*\[Tau]*theta]; Mat[[i,i + M]] = -h^2; ]
For[i = M + 1, i <= 2*M, i++, If[i != M + 1, Mat[[i,i - 1]] = (-beta)*\[Tau]*theta]; Mat[[i,i]] = h^2 + 2*beta*\[Tau]*theta; 
   If[i != 2*M, Mat[[i,i + 1]] = (-beta)*\[Tau]*theta]; Mat[[i,i - M]] = -h^2; ]

solver = LinearSolve[Mat]; 

For[i = 0, i <= M + 1, i++, u[i, 0] = u0[i*h]; v[i, 0] = v0[i*h]]

For[n = 1, n <= Nd, n++, 
  For[i = 1, i <= M, i++, ps[[i]] = h^2*(u[i, n - 1] - v[i, n - 1]) + (1 - theta)*\[Tau]*alpha*(u[i - 1, n - 1] - 2*u[i, n - 1] + 
         u[i + 1, n - 1]); ps[[M + i]] = h^2*(v[i, n - 1] - u[i, n - 1]) + (1 - theta)*\[Tau]*beta*(v[i - 1, n - 1] - 2*v[i, n - 1] + 
         v[i + 1, n - 1])]; r = solver[ps]; u[0, n] = u[M + 1, n] = v[0, n] = v[M + 1, n] = 0.; 
   For[i = 1, i <= M, i++, u[i, n] = r[[i]]; v[i, n] = r[[M + i]]]; ]

n = 0; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]
n = 1; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]
n = 5; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]

其它你可能感兴趣的问题

上一篇信任域算法的全局收敛下一篇缺少初始条件的 PDE 类有名称吗？