2

我正在尝试解决具有不确定形式的微分方程,我正在使用 Julia 和 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)

这个问题的解析解(在天体物理学中被称为邦迪流问题)是众所周知的(即使对于不确定的情况)。我对从求解器获得的不确定解感兴趣。

当我使用求解时,sol = solve(prob)我得到的振荡解与我知道的平滑解析解完全不同(参见下面链接中的图)。

绘制 u 与 t

我期待在 t 接近 50 时遇到一些“问题”(同时,由 u 表示的 y 轴变量(表示速度)将接近 100),因为只有这样分子(和分母)才会一起消失。任何想法为什么解决方案开始振荡?

我也尝试过sol = solve(prob, alg_hints = [:stiff])并收到以下警告:

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

解决方案仍然振荡(类似于在没有施加刚度的情况下获得的解决方案)。

我在这里做错了什么吗?是否有另一种方法可以使用DifferentialEquations.jl 包解决此类不确定方程?

4

1 回答 1

0

该 ODE 在奇点处在数值上非常不稳定。遇到此类问题时,基于时间步进的方法很容易爆炸。您需要使用基于搭配的方法,例如 ApproxFun.jl 中的某些方法,以避免直接评估奇点以稳定求解。

于 2021-05-21T12:25:51.897 回答