扩散方程的C-N格式计算问题

hxl166417

大家好
        一般的扩散方程例如：∂u/∂t=a·∂*2u/∂x*2  （一维） ，当常系数a的值为+值时，我编写的C-N格式求解没有问题。但是，当系数a的值为负值时，用原本的C-N格式求解会出现数值震荡，我知道系数为负值没有物理意义，只是想求解看看。
      请问有没高手知道原因。多谢！

onesupeng

是发散还是数值震荡？这可是一个反扩散的过程，很快就无穷大了阿

hxl166417

计算几步之后，数值就无穷大了，应该叫“数值色散”吧？反扩散过程是一个什么样的物理过程？我用的初始化条件是个“矩形波”。在数值突变的地方，数值震荡越来越大。

   谢谢您的指点。

lwd1981

这不是一个困难的问题。你的耗散是负的，意味着系统的能量随着时间是增加的，最终必然无穷大。
建议你做一下几件事情：
1.你可以用分离变量方法求一个解析解看看
   U=f(t)*g(x)来求解析解，看看是不是随着时间增大，很快变得无穷大；
2.  你也可以直接将上述方程两边乘以u，然后积分变成左侧是能量的2范数随时间的变化，右侧是一个东西（加上边界条件），你看看系数a对系统能量的影响。
3. 你还可以将方程数值离散后，无论采用能量法，还是傅立叶分析方法看你稳定性对系数a的要求；

hxl166417

谢谢您的指点

我用分离变量法求出了解析解。U数值随着时间的增大，确实会变得无穷大。但是这跟计算过程中的数值震荡有联系吗？我用的初始化条件是个“矩形波”。在数值突变的地方，数值会出现震荡。

2和3里面的知识我不是很懂。我再学习一下，再次感谢您的帮助。

lwd1981

原帖由 hxl166417 于 2014-2-21 08:02 发表

                               
登录/注册后可看大图

谢谢您的指点

我用分离变量法求出了解析解。U数值随着时间的增大，确实会变得无穷大。但是这跟计算过程中的数值震荡有联系吗？我用的初始化条件是个“矩形波”。在数值突变的地方，数值会出现震荡。

2和3里面的 ...

你可以将你的数值离散格式对应的修正方程导出来，你会发现你的数值耗散项也是负的；如果你的数值耗散是正的，那么出现间断附近（数值突变的地方）的振荡将被光滑掉（这其实就是“人工粘性”的机制），从而使得你的“矩形波”会变光滑；但是你的数值耗散项是负的，它不会将会将间断附近（数值突变的地方）所出现的振荡光滑，反而会将间断附近（数值突变的地方）所出现的振荡放大。

至于2实际上也是能量法。根据能量最小原理，一个稳定系统的能量随着时间的变化率应该小于或者等于0. 能量法是一个非常好的用于判断系统是否稳定的方法。一般而言，目前有关偏微分方程数值格式的稳定性分析大都是针对线性方程或者拟线性方程进行分析，直接对于非线性方程目前可用的方法不多，能量法就是一种。如果你对高阶方法感兴趣，你会发现能量法几乎成为分析高阶方法稳定性的唯一方法了。实际上，能量法不仅可以用来分析稳定性，还可以分析应该如何让施加边界条件。《nodal discontinuous galerkin methods:algorithms analysis and applications》一书中对能量法有详细的应用阐述。

至于3是数值分析里面用来分析一个数值格式的稳定性以及稳定区的基本方法，建议你看看有关偏微分方程数值计算方面的书籍。

[ 本帖最后由 lwd1981 于 2014-2-21 12:41 编辑 ]

coolboy

原帖由 lwd1981 于 2014-2-21 20:28 发表

                               
登录/注册后可看大图

你可以将你的数值离散格式对应的修正方程导出来，你会发现你的数值耗散项也是负的；如果你的数值耗散是正的，那么出现间断附近（数值突变的地方）的振荡将被光滑掉（这其实就是“人工粘性”的机制），从而使得你的“矩形波”会变光滑；但是你的数值耗散项是负的，它不会将会将间断附近（数值突变的地方）所出现的振荡光滑，反而会将间断附近（数值突变的地方）所出现的振荡放大。

这一段的叙述表明发言者对问题的数学物理内涵有着比较清晰的认识，即数值格式通常在解的“间断”处会产生虚假波以及稳定的数值格式会出现“人工粘性”来耗散掉（物理的或虚假的）数值波动。

现在我来提出一个问题，看看是否有人知道答案或知道的是什么样的答案。假如实际的物理问题确实存在类似于上述“矩形波”即具有间断的特性但方程中并不存在物理的“粘性耗散”过程，则该如何设计数值格式使得守恒的“矩形波”特性在经过长时间或长距离积分之后仍然能被很好地保持着。

至于实际的物理过程，例如下帖子中提到的守恒的总能量若初始时刻刚好是一个矩形波，则它就可以被实际流场传输很远或很久而大致不变形，至少理论上是如此。现在的问题是数值求解此问题时该如何实现？

+++++++++++++++
[讨论]伯努利方程是能量方程还是动量方程？
http://www.cfluid.com/bbs/viewthread.php?tid=114265
+++++++++++++++

数学家冯康后期有一出名的工作就是数值求解哈密顿系统并保持其保守性。但他处理的问题是常微分方程（组）而我们现在的问题是要数值求解双曲型的偏微分方程（组）并使得解的“数值粘性耗散”尽可能地小。

hxl166417

您说的很有道理，让我清晰明了的认识到这个问题，非常感谢！我能不能这样理解，现在出现数值震荡问题，是这个方程本身的性质决定的。不能改变了。
      我还有两个问题：1.这个方程最后的计算结果肯定是无穷大，但是在数值计算过程中就出现了数值震荡，是不是正常？2.我要是添加“人工粘性项”来消除这个震荡，对最终的结果的合理性会有“本质影响”吗？
     您的回到让我感受到自己需要补充的知识还很多，多谢您的指点！

lwd1981

原帖由 hxl166417 于 2014-2-22 08:35 发表

                               
登录/注册后可看大图

您说的很有道理，让我清晰明了的认识到这个问题，非常感谢！我能不能这样理解，现在出现数值震荡问题，是这个方程本身的性质决定的。不能改变了。
      我还有两个问题：1.这个方程最后的计算结果肯定是无穷大，但 ...

首先你这里的震荡不是方程本身的性质决定，引起震荡是由于你给出的方波形初值造成的。我只是说你所采用的计算方法不能消除这个震荡罢了！
   你可以通过调节人工粘性项来消除震荡，这没有问题，关键问题在于添加多大的人工耗散！这个多大不好拿捏！需要设计精巧的数值耗散，一不影响数值精度，二可以有效消除震荡。

通流

我有点糊涂了。初值怎么是造成震荡的原因？应该是这个格式对于一个间断造成的震荡。也就是问题是格式造成的。不过好像大部分格式都有这个问题，所以才需要一点点数值粘性。

glandetian

空间上采用WENO重构，时间上用辛算法？
据报道，辛算法在处理长时间问题有优势，如在卫星编队飞行、星际航行等方面有应用。据钟万勰院士的意思，辛算法在保辛的同时，也可以保能量。以前的看法是保辛就不保能量，保能量就不保辛。
如果方程是线性的，在时间上可以用精细积分方法。这个算法可以确保长时间模拟的波形保真，因为这个算法在时间方向几乎就是解析解。当然，计算规模不能太大，因为精细积分涉及到20次矩阵与矩阵的相乘计算，规模大的话，内存和计算时间上都受不了。如果方程是非线性的，精细积分可能就不行了，虽然有人提出了非线性常微分方程的増维精细积分方法，但是，毕竟已经不能接近“解析解”，所以，长时间模拟就成问题了。
当然，WENO重构本来就是非线性的，所以。。。。呵呵，白说了

coolboy

原帖由 glandetian 于 2014-2-23 22:40 发表

                               
登录/注册后可看大图

空间上采用WENO重构，时间上用辛算法？
据报道，辛算法在处理长时间问题有优势，如在卫星编队飞行、星际航行等方面有应用。据钟万勰院士的意思，辛算法在保辛的同时，也可以保能量。以前的看法是保辛就不保能量，保 ...

我不知道什么是“WENO重构”。辛算子是指算子的空间结构，“时间上用辛算法”是什么意思？

钟万勰的工作是在结构力学和控制论领域，故他的算法也应该是关于常微分方程组而非流体力学的偏微分方程。

据钟万勰的意思，哈密顿系统或辛算子也就是中国古代的“八卦两仪”或与之密切联系。我对此持保留观点或认为他是在瞎扯，呵呵。据此（或即依此类推、或根据他对科学研究的不严肃态度），我对他那所谓“辛算法在保辛的同时，也可以保能量”什么的意思也持保留态度，呵呵。

glandetian

我最开始看到“辛算法”，是在计算电磁学上，如吴先良基于辛算法的二维电磁场散射问题的研究。后来，在余德浩的《微分方程数值解法》里看到有辛Runge-Kutta方法的介绍，可以用辛算法来做时间方向上的离散，保持长期追踪计算的准确性。并没有更深入的学习，可能有说错的地方，呵呵。。。见谅啊！ 不过，从吴先良的论文看，确实空间离散与时间离散是分开的，然后，时间离散用辛格式。
WENO是“加权本质无振荡”的缩写，是在做激波捕捉时提出来的，计算效果很好。在用差分法计算Euler方程时，WENO格式直接包含在计算过程中；在用间断Galerkin方法求解Euler方程时，在每个时间步结束后（或者在每个时间步的每个分步里面）做WENO重构。我想，您说的矩形波，应该跟激波差不多意思吧，所以想到了WENO格式。
流体力学方程组，在空间离散完后，也得到一个非线性常微分方程组啊，所以，可以借鉴钟万勰提出的关于常微分方程组的数值解法。
对于钟万勰的“保辛同时保能量”说法，我还不太懂，呵呵。。。

静候您对自己提出的问题的答案！  ^_^

[ 本帖最后由 glandetian 于 2014-2-24 07:44 编辑 ]

onesupeng

好奇地问一下，谁能写出这个问题的解析解？

lwd1981

原帖由 onesupeng 于 2014-2-23 23:59 发表

                               
登录/注册后可看大图

好奇地问一下，谁能写出这个问题的解析解？

由于不知道他具体的边界条件，不过对于1为纯扩散或者热传导，只要有就界就可以求出解析解。下面就是一个简单边界条件的解法：
http://zh.wikipedia.org/wiki/%E7%86%B1%E5%82%B3%E5%B0%8E%E6%96%B9%E7%A8%8B%E5%BC%8F