如何用DifferentialEquations.jl消除不定形式微分方程解的振荡

Question

我试图用不确定的形式来解微分方程，我正在玩朱莉娅和DifferentialEquations.jl软件包来解决这个问题。

我从一个不确定形式的简单微分方程开始(作为一个测试用例)，看看它是否可行(我的最终目标是解决一个复杂的微分方程组，并将它与机器学习结合起来。这需要我进行大量的分析和编码，我只想在我知道一些基本的测试用例起作用之后才开始)。

下面是我开始使用的测试用例代码：

bondi(u,p,t) = -(2*k/t^2)*(1 - (2*a^2*t/k))/(1 - a^2/u)
u0 = 0.01
p = (1e4, 10)
k,a = p
tc = k/(2*a^2)
tspan = (tc/10, tc*5)
prob = ODEProblem(bondi,u0,tspan,p)

这个问题的解析解(在天体物理学中称为Bondi流问题)是众所周知的(即使在不确定的情况下)。我感兴趣的是从求解者那里得到的不定解。

当我使用sol = solve(prob)进行求解时，我得到了一个与我所知道的光滑解析解完全不同的振荡解(参见下面链接中的图)。

图u对t

当t接近50时，我预期会遇到一些‘问题’(同时，u表示的y轴变量(表示速度)接近100)，因为只有这样分子(和分母)才会一起消失。知道为什么这些解决方案会开始振荡吗？

我还尝试过使用sol = solve(prob, alg_hints = [:stiff])，并得到了以下警告：

警告:中断。需要更大的上限。@ DiffEqBase C:\Users\User.julia\packages\DiffEqBase\1yTcS\src\integrator_interface.jl:329

解仍在振荡(类似于在不施加刚度的情况下得到的解)。

我在这里做错什么了吗？还有另一种方法可以用DifferentialEquations.jl软件包来解这种不定方程吗？

Answer 1

在奇异点，该ODE在数值上非常不稳定。基于时间步长的方法在遇到这种问题时很容易爆炸。您需要使用一种基于配置的方法，如ApproxFun.jl中的某些方法，以避免在奇点处直接进行评估，以便稳定求解。