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

Question

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

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

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

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

Спасибо заранее

Francisco

Answer 1

Ответ на кросс-пост ОП на Julia Discourse ; скопировано здесь для полноты.

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