0

我是 Julia 的新手，我想解决这个系统：

在此处输入图像描述

其中 k1 和 k2 是常数参数。然而，当 y,0 时 I=0，否则为 Ky，其中 k 是一个常数值。

我按照关于 ODE 的教程进行操作。问题是，如何在DifferentialEquations.jl 中求解这个分段微分方程？

4

1 回答 1

1

在Julia Discourse上的 OP 交叉帖子上回答；为了完整起见，复制到这里。

这是一个（稍微）有趣的例子 $x''+x'+x=\pm p_1$，当在 $x=p_2$ 处遇到切换流形时，$p_1$ 的符号会发生变化。为了让事情变得更有趣，请考虑开关歧管中的滞后，这样每当穿过开关歧管时，$p_2\mapsto -p_2$。

代码相对简单；StaticArrays/SVector/MVector 可以忽略，它们只是为了速度。

using OrdinaryDiffEq
using StaticArrays

f(x, p, t) = SVector(x[2], -x[2]-x[1]+p[1])  # x'' + x' + x = ±p₁
h(u, t, integrator) = u[1]-integrator.p[2]  # switching surface x = ±p₂;
g(integrator) = (integrator.p .= -integrator.p)  # impact map (p₁, p₂) = -(p₁, p₂)

prob = ODEProblem(f,  # RHS
                  SVector(0.0, 1.0),  # initial value
                  (0.0, 100.0),  # time interval
                  MVector(1.0, 1.0))  # parameters
cb = ContinuousCallback(h, g)
sol = solve(prob, Vern6(), callback=cb, dtmax=0.1)

然后绘图sol[2,:]以sol[1,:]查看相平面 - 在这种情况下是一个很好的非光滑极限环。

请注意，如果您尝试使用结果解（即sol(t)）的插值，则需要非常小心具有不连续导数的点，因为插值有点错误。这就是为什么我习惯dtmax=0.1在这种情况下获得更平滑的解决方案输出。（我可能也没有使用最合适的积分器，但它是我在之前复制和粘贴的代码中使用的积分器。）

于 2018-07-05T23:33:50.223 回答