如何使用MATLAB求解这个偏微分方程

Question

我在一份实践考试中遇到了以下问题：

我需要使用MATLAB进行求解。问题是，我以前没有遇到过这样的问题，很难开始解决。

我有一个1x1的网格，分成10x10。我知道我可以使用1/10*x*2计算除角落以外的整个底行。我还知道我可以使用(1/10)(1+t)^2计算整个右侧行。但是，我无法弄清楚如何获得足够的点来填写整个网格的值。我知道这一定与问题中给出的偏导数有关，但我不太确定它们如何发挥作用（尤其是u_x方程）。有人能帮我入门吗？

我不需要整个解决方案。一旦我有足够的点，我就可以轻松编写Matlab程序来解决其他问题。实际上，我认为我只需要解决x=0轴，然后填写网格的中间部分即可。

我已经计算出了底行（两个角除外），分别为0.001，0.004，0.009，0.016，0.025，0.036，0.049，0.064，0.081。同样地，整个右边行根据给定的边界条件很容易计算。我只是不能将它们拼凑在一起。

编辑：第三个边界条件方程式拼写错误。它应该是： u_x(0,t) = 1/5t，而不是u(0,t) = 1/5t。

- user3330644

Answer 1

首先要认识到你需要解决的方程是“线性波动方程”，而你所获得的数值方案可以重写为：

( u^(n+1)_m - 2u^n_m + u^(n-1)_m )/k^2 = ( u^n_(m-1) - 2u^n_m + u^n_(m+1) )/h^2

其中k是时间步长，h是空间中的delta x。

重新构造的数值方案表明左侧和右侧是二阶中心有限差分逼近u_tt和u_xx。

然而，为了在数值上解决问题，您需要使用给定的形式，因为它是您需要在数值上实现的显式更新公式：它将时间n+1的解作为前两个时间n和n-1的函数给出。您需要从初始条件开始，并按时间推进解决方案。

请注意，解决方案在域的边界（x=0和x=1）上被分配，因此离散化解决方案u^(n)_0和u^(n)_10的值对于任何n（t=n*k）都是已知的。在第n个时间步骤中，您的未知向量是[u^(n+1)_1, u^(n+1)_2, ..., u^(n+1)_9]。

请注意，使用更新公式来找到第n+1步的解，需要知道前两步的解。那么，如果您需要从前两个时间获取信息，如何从n=0开始？这就是初始条件发挥作用的地方。

您在n=0（t=0）处有解，但您还在t=0处有u_t。这两个信息结合起来可以给您提供u^0和u^1并让您开始。

我将使用以下启动方案：

u^0_m = u(h*m,0)  // initial condition on u
(u^2_m - u^0_m)/(2k) = u_t(h*m,0)  // initial condition on u_t

结合使用 n=1 的数值方案，您可以定义 m=1,...,9 的 u^1_m 和 u^2_m 的线性系统所需的一切。

总之：

--使用启动方案同时找到 n=1 和 n=2 的解。

--从那里开始使用给定的数值方案推进时间。

如果您完全迷失，请查看有关：有限差分方案、对流方程的有限差分方案、双曲型方程的有限差分方案、时间推进等内容。

编辑：

对于 u_x 的边界条件，通常使用 ghost cell 方法：

• 在m=-1处引入幽灵单元，即虚构（或辅助）网格点，用于处理边界条件，但不属于解决方案。

• 第一个节点m=0回到未知向量中，即现在使用的是[u_0 u_1 ... u_9]。

• 使用左侧边界条件来关闭系统。具体而言，通过编写边界条件的中心近似式

  u^n_(1) - u^n_(-1) = 2*h*u_x(0,k*n)

• 上述方程允许您以内部真实节点的解的形式表达幽灵节点上的解。因此，您可以将时间推进数值方案（给定的方案）应用于m=0节点。（应用于m=0的数值方案将包含来自m=-1幽灵节点的贡献，但现在您已经将其表示为m=1节点的形式。）

- Emerald Weapon

---

我犯了一个只有一个字符的错别字，可能会稍微改变答案。由于我不知道如何自己输入这些方程式，所以我进行了编辑来修复它。第三个边界条件不是直接赋值的边界条件。这个边界条件让我很困扰，因为它阻止了我获得左列的值，而我需要用目前已有的知识来解决这个问题。 - user3330644

请查看我回答的相应编辑。 - Emerald Weapon

所以，我得到了n=0行和m=1列。然后我使用u_t同时计算n=1和n=2？我不确定那是怎么回事。由于u_t方程只指定了t为0时，我能不能只用它来找到n+1层呢？ - user3330644