0

I am new to Julia, I would like to solve this system:

enter image description here

where k1 and k2 are constant parameters. However, I=0 when y,0 or Ky otherwise, where k is a constant value.

I followed the tutorial about ODE. The question is, how to solve this piecewise differential equation in DifferentialEquations.jl?

JohanC
  • 71,591
  • 8
  • 33
  • 66
Francisco
  • 1

1 Answers1

1

Answered on the OP's cross post on Julia Discourse; copied here for completeness.

Here is a (mildly) interesting example $x''+x'+x=\pm p_1$ where the sign of $p_1$ changes when a switching manifold is encountered at $x=p_2$. To make things more interesting, consider hysteresis in the switching manifold such that $p_2\mapsto -p_2$ whenever the switching manifold is crossed.

The code is relatively straightforward; the StaticArrays/SVector/MVector can be ignored, they are only for speed.

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)

Then plot sol[2,:] against sol[1,:] to see the phase plane - a nice non-smooth limit cycle in this case.

Note that if you try to use interpolation of the resulting solution (i.e., sol(t)) you need to be very careful around the points that have a discontinuous derivative as the interpolant goes a little awry. That's why I've used dtmax=0.1 to get a smoother solution output in this case. (I'm probably not using the most appropriate integrator either but it's the one that I was using in a previous piece of code that I copied-and-pasted.)

DavidB
  • 36
  • 1