计算科学 - 在 matlab 中使用 ODE 选择步长 - 吾爱随笔录

在 matlab 中使用 ODE 选择步长

计算科学 matlab 颂量子力学

2021-11-26 02:52:25

嘿，谢谢你抽出时间来看看我的问题。这是我之前在physics.stackexchange.com 上发布的问题的更新版本

我目前正在研究 2D 激子旋量 Bose-Einstein Condensate，并对这个系统的基态感到好奇。到达基态的数学方法称为虚时间法。

该方法非常简单，将量子力学中的时间替换为虚数这种替换导致我系统中的高能粒子比低能粒子衰减得更快。在计算的每一步重新归一化粒子数量，我们最终得到了一个能量最低的粒子系统，也就是。基态。

t = - i τ

$t = -i \tau$

所讨论的方程是非线性的，称为非线性薛定谔方程，有时称为Gross-Pitaevskii 方程。为了解决这个问题，我使用了 Matlabs ode45，它使系统及时向前发展并最终达到基态。

笔记！非线性薛定谔方程包括空间中的拉普拉斯和其他一些微分项。这些都是使用快速傅里叶变换解决的。最后我们只有一个时间 ODE。*

我的问题和疑问：计算从到。ode45 被放入一个for循环中，因此它不会同时计算一个巨大的向量。第一轮将从 ode45(odefun, ) 开始，然后下一次从开始。这里时间步长是我的问题。时间步长的不同选择给了我不同的基态解决方案，我不知道如何确定哪个时间步长给了我“最”正确的基态！ $t_0$ $t_f$ $[t_0,\dots,t_f]$ $[t_0, t_0+\Delta/2, t_0 + \Delta],y,\dots$ $t_0 + \Delta$ $\Delta$

我的尝试：我意识到在这个方案中，大时间步长会导致大量粒子在重新归一化为原始粒子数之前衰减，而小时间步长会导致更少量的粒子在重新归一化之前衰减。我最初的想法是小时间步长应该提供更准确的解决方案，但似乎恰恰相反。

我不是数字专家，所以选择 ode45 只是随意的。ode113 给了我同样的东西。:(

有没有人对这个问题有任何想法。让我知道是否需要任何额外的细节。

谢谢你。

更新 1： 我一直在研究虚时间方法和 ODE。似乎如果时间步长不够小，整个事情就会变得不稳定。这让我想知道我的非线性方程是否僵硬，这使我理解的事情变得更加困难。我会及时通知你的。

更新 2： 已修复：问题确实是在 ODE 之外进行了规范化。如果归一化保留在 odefun 中，则 ODE 针对“外部”时间步的不同选择返回相同的结果。我的同事向我展示了旧代码，我只是在我的 odefun 中添加了一行。

function y_out = odefun(t,y_in,...variables...) 

    ...
    [ Nonlinear equations evaluated ]  
    ...


    y_out = y_out + 0.1*y_in*(N0-Ntemp) ;
end

最后一行计算当前粒子数 (Ntemp) 和系统应容纳的粒子数 (N0) 的差值。它将一部分粒子添加回输出，从而在系统中创建总粒子数稳定性，而不是让它们全部衰减。

我还将提出一个关于问题维度的新问题，以及在 ODE 中使用皮秒或纳秒作为时间步长的一些差异。

谢谢你们。:)

3个回答

由于您没有发布 MATLAB 代码，我不确定您是如何调用 ode45 的。我猜您在每次调用 ode45 时都会更改 tspan 向量（第二个参数）。首先要了解的是 tspan 向量（几乎）对 ode45 使用的时间步长没有影响。tspan 向量只允许您将积分的时间跨度以及您想要输出的时间传递给 ode45。ode45 中 Runga-Kutta 算法使用的时间步长在内部进行了调整，以达到规定的精度。控制此精度的两个参数是传递给 ode45 的选项结构中的 RelTol 和 AbsTol。它们具有合理的默认值，并且由于您没有提及这些，我假设您没有更改它们。

我说“几乎”对正常的 ode45 时间步没有影响。如果您以非常小的时间间隔请求输出，相对于 ode45 否则将采取的时间步长，那么它将必须减少时间步长以满足您的输出请求。我相信这就是 JM 假设正在发生的事情。为什么你需要在这么多输出时间的解决方案？通常只需在足够多的时间请求输出以生成平滑图就足够了。

至于您看到的解决方案的变化，可能 RelTol 和 AbsTol 的默认值不适合您的问题。我建议用一次调用替换 ode45 上的循环，以合理的次数请求输出，并尝试使用较小的 RelTol 和 AbsTol 值，直到获得收敛的解决方案。

由于非线性薛定谔方程是非线性的，它可以有大量的稳态，其中一些可能是稳定的。在物理现实中，从一种特定状态开始，系统将确定性地演变为一种最终状态。如果数值方案为不同的离散化（时间步长）提供不同的结果，那么这是离散化的根本缺陷。检查你的代码。

如果你最终得到一个状态，那么很容易验证它是否真的是一个静止状态：如果时间演化由然后如果你确实最终得到不同的静止状态，你可以比较它们的吉布斯能量其中是能量密度。当时，通常看起来很简单，例如，对于 Ginzburg-Landau 方程，。 $\psi_0$

\frac{d ψ}{d t} = F (ψ),

$\frac{\text{d}\psi}{\text{d}t} = F(\psi),$

F (ψ_{0}) = 0.

$F(\psi_0) = 0.$

G (ψ) = \int_{Ω} E (ψ)

$G(\psi) = \int_\Omega E(\psi)$

E (\cdot)

$E(\cdot)$

F (ψ) = 0

$F(\psi)=0$

E (ψ)

$E(\psi)$

E (ψ) = - | ψ |^{4}

$E(\psi) = -|\psi|^4$

问题解决了：

规范化需要成为 ODE 中评估的函数的一部分。在多个步骤中分解 ODE 并在它们之间进行归一化会导致看似数值不稳定，并根据 ODE 分解的时间间隔产生不同的结果。（有关更多详细信息，请参阅相关编辑 2。）

其它你可能感兴趣的问题

上一篇八度：计算两个向量矩阵之间的距离下一篇当解处的雅可比行列式为奇异时牛顿法的策略