Кусочно-дифференциальное уравнение с использованием Юлии

Question

Я новичок в Юлии, я хотел бы решить эту систему

\frac{dx}{dt} = k1y \ \frac{dy}{dt} = k2y+I

где k1 и k2 - постоянные параметры. однако, I=0, когда y,0 или Ky в противном случае, где k является постоянным значением.

Я следовал учебному пособию по ODE в вопросе, как можно решить это кусочно-дифференциальное уравнение в DifferentialEquations.jl, однако я не нашел ничего о кусочно-зависимых функциях

Заранее спасибо

Франциско

Answer 1

Ответил на перекрестный пост ОП на Джулии Дискурс; здесь скопировано для полноты.

Вот (слегка) интересный пример $x''+x'+x=\pm p_1$, где знак $p_1$ изменяется, когда встречающийся многообразие встречается при $x=p_2$. Чтобы сделать вещи более интересными, рассмотрим гистерезис в переключающем многообразии, так что $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 чтобы получить более плавный вывод решения в этом случае. (Возможно, я тоже не использую наиболее подходящий интегратор, но это тот, который я использовал в предыдущем фрагменте кода, который я скопировал и вставил.)