在自适应步长 SDE 积分器中测试 Wiener 过程拆分

计算科学 测试 随机抽样 随机颂歌 自适应时间步长
2021-11-26 08:07:34

我正在研究用于随机微分方程的自适应步长积分的各种方法并尝试实现它们。我看过的所有论文（例如H. Lamba, J. Comp. App. Math. 161(2), 417–430 (2003)）都以某种形式提到必须始终保留维纳过程。换句话说，如果某些步骤被拒绝，仍然必须检查所有稍后定位的采样值。

例如，假设你在时间t0并试图走到t1，采样维纳增量dW1为了那个原因。该步骤结果出现不可接受的错误并被拒绝。现在你尝试迈出一步t2<t1. 看来，正确的程序是对中间增量进行采样dW2从t0到t2基于dW1（确切的公式在这里并不重要），并从t2到t1成为dW1−dW2.

现在假设这一步也被拒绝了。重复上面的过程，你现在有三个增量，dW3,dW2−dW3和dW1−dW2从步骤t0到t3，然后到t2然后到t1.

这部分我没有看到解释，虽然似乎所有作者都暗示如果你的步骤t3成功了，你的下一步一定不能过去t2并且应该使用采样增量dW2对于那部分。如果你遵循这个逻辑，在积分器的实现中，你将不得不维护一个采样增量的动态列表，这会增加内存占用（也以不可预测的方式）。据推测，这避免了采样维纳过程中的偏差。我们称之为方法 A。

我试图用更简单的逻辑替换它，其中只有最远的时间样本（dW1) 被保留。这样，在上述情况下，您采样后dW3， 你忘记dW2，并且对于您的下一步，采样增量仅取决于dW1. 我们称之为方法 B。

抱歉，解释太长了，现在回答我的问题。假设我的简化逻辑应该在解决方案中引入偏见。问题是我无法在我尝试的任何测试中检测到它。您对如何设计能够区分方法 A 和方法 B 的测试有任何想法吗？

1个回答

我将在本文（Rackauckas 和 Nie 2017）中更详细地讨论您描述为 RSwM2 的方法。在那篇论文中，我稍微能够检测到它有时会做错什么，但由于它只有被拒绝的问题，所以没什么大不了的。这 3 种方法（RSwM1、RSwM2、RSwM3）现在是DiffEqNoiseProcess.jl的基础，随后是在DifferentialEquations.jl中构建自适应时间步长的方法。

不用说，我现在的测试比论文发表时要好得多。首先要注意的是，如果您切换到使用 PI 控制的自适应时间步长（详细信息请参阅 Hairer II），这将不会成为问题，因为这会极大地减少重新拒绝的机会。调整自适应参数也有帮助。但是，一个很好的测试是使用自适应时间步来概括布朗桥的特性。使用论文中的符号，如果您从将值放入未来信息堆栈开始S1在最后的时间点，你必须达到那个值。所以这意味着从一个值开始S1应该使噪声过程成为布朗桥。由于您知道布朗桥在每个时间点的均值和方差等属性，因此很容易运行大量蒙特卡罗模拟并以非常敏感的方式测试这些属性。我在我的持续集成测试套件中有这个测试，并且发现它对让算法正确非常敏感。为了检查，我使用 RSwM2 实现运行它：

W = BrownianBridge(0.0,1.0,0.0,1.0,0.0,0.0;rswm = RSWM(adaptivealg=:RSwM2))
prob = NoiseProblem(W,(0.0,1.0))
monte_prob = MonteCarloProblem(prob)
@time sol = solve(monte_prob,dt=0.1,num_monte=10000)

# Spot check the mean and the variance
qs = 0:0.1:1
for i in 2:10
  q = qs[i]
  @test ≈(timestep_mean(sol,i),q,atol=1e-2)
  @test ≈(timestep_meanvar(sol,i)[2],(1-q)*q,atol=1e-2)
end

这个测试仍然是随机的，但通过 1 和 3 的次数往往多于 2。它的容差是这样设置的，这样它就可以在 Travis 和 Appveyor 上运行得足够快，但如果你增加数量，它应该能够更清晰地区分它们轨迹。

当然，您可以在自己的实现中做同样的事情来查看。如果你能得到一个具有突然刚度以可靠地诱导再拒绝的解析解，那么测试可能会变得更好，但至少布朗桥变体是我迄今为止最严格的测试。

其它你可能感兴趣的问题
上一篇GMSH 中的自定义网格尺寸 下一篇解决具有相同约束的多个最小二乘问题