algorithm - 如何使用 MATLAB 求解此 PDE

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

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您未知的第 th 时间步是矢量 [u^(n+1)_1, u^(n+1)_2, ..., u^(n+1)_9] .

还要注意使用更新公式在 n+1 处找到解决方案步骤，需要了解前面两个 步的解决方案。那么，你如何从n=0开始呢？如果您需要以前两次 的信息？这是初始条件发挥作用的地方。您可以在 n=0 找到解决方案( t=0 )，但你也有 u_t在t=0 .这两条信息结合起来可以给你两个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 使用的数字方案相结合为您提供为 u^1_m 定义线性系统所需的一切和 u^2_m对于 m=1,...,9 .

总结:

--使用启动方案在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 节点表示。)