0

I'm teaching Differential Equations and am using Pluto notebooks with DifferentialEquations.jl to supplement the course. We are currently trying to model a spring-mass system that is subject to an impulse of 10 Newton-seconds (hammer blow) at a specific instant of time. I believe I need to use a callback to simulate this, but don't know how to implement it.

The differential equation for this system is: u''(t) + 2u'(t) + 10u(t) = 10*delta(t-1), where 10*delta(t-1) is the delta function for the 10 Newton-second impulse at 1 second. The initial conditions are u(0) = u'(0) = 0.

What I've done so far is:

function hammer_time(du, u, p, t)
    -2*du - 10*u
end

begin
    du0 = 0.0
    u0 = 0.0
    tspan = (0.0, 10.0)
end

prob = SecondOrderODEProblem(hammer_time, du0, u0, tspan)

begin
    hit_times = [1.0]
    affect!(integrator) = integrator.u += 0.5

    cb = PresetTimeCallback(hit_times, affect!)
    sol = solve(prob, Tsit5(), callback = cb)
end

In the last code block I know the second line is wrong, but I don't know what it should be to model the hammer blow impulse. I assume it changes the energy state of the system but I don't know what to do to model that.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
G. Church
  • 15
  • 3
  • Please read again the interface description. `du` is the output of the ODE function, `u` the input state. Usually `u[0]` stands for the solution value, `u[1]` for the first derivative, and the output value `du[0]` then should be the second derivative. But compare this idea again against the documentation of all parts. // The delta transformed as action should cause a jump by 10 in `u[1]`, the derivative value. // The first part is insufficient study of the documentation and thus borderline off-topic, the second part is borderline on-topic in what concerns the correct use of timed actions. – Lutz Lehmann Apr 29 '23 at 17:28
  • I appreciate the reply, but, frankly, I don't find it very helpful. I was hoping someone could show me how to modify the code sample I provided, to get the effect I want. I have looked through [the documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/), but being inexperienced with callbacks and not a professional programmer, I find it somewhat opaque. – G. Church Apr 30 '23 at 02:35
  • 1
    It was not meant to be directly helpful. It was a prompt to correct the absolute basics, to ensure you start from a minimal example without actions/affectors that compiles and works. And also to review the mathematical side of it, how to translate non-representable distributions into implementable effects. Revise the guidelines for well-posed questions. – Lutz Lehmann Apr 30 '23 at 07:53

1 Answers1

2

As far as I can tell, SecondOrderODEProblem does not play nicely with callbacks mutating states - it originates from some details under the hood when constructing this nice interface, not important for your use case though.

However, we can certainly avoid this altogether if we write the 2nd order ODE as a system of first order ODEs:

Let u = [x; v], so that u' = [v; -2v-10x]

Then we can define our system using the more common ODEProblem interface

using DifferentialEquations, Plots;

#x = u[1], v = u[2]
function hammer_time(u, p, t)
    return [u[2];
            -2*u[2] - 10*u[1]];
end

u0 = [0.0;
      0.0];
tspan = (0.0, 10.0);

condition(u, t, integrator) = (t == 1);
affect!(integrator) = (integrator.u[2] += 10);
cb = DiscreteCallback(condition, affect!);
prob = ODEProblem(hammer_time, u0, tspan)
sol = solve(prob, Tsit5(), callback=cb, tstops=[1.0])
plot(sol, labels=["x" "v"])

This has the effect that at time 1, we increase the velocity (u[2]) by 10 units.

The above code yields plot of position and velocity as a function of time.

I'm sure you already found the documentation for callbacks, but here it is again for posterity. https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#Using-Callbacks

Adrian Mole
  • 49,934
  • 160
  • 51
  • 83
Criston
  • 36
  • 4