In Julia, I want to use DifferentialEquations.jl package to solve
\ddot{u} + f(u,\dot{u},p) = g(t)
where g(t) is in given as a vector of values at equidistanced instants of time t.
This situation is different from one in https://diffeq.sciml.ai/stable/tutorials/ode_example/ where the forcing function M(t) is continuous.
A solution to this situation is to interpolate g(t) as suggested in
solve system of ODEs with read in external forcing
However, I do not want to interpolate the forcing function g(t) but want to try callback command instead.
For a 2-order linear SDOF system subjected to ground motion, I already tried
using DifferentialEquations
rhs = rand(10); # ground accel.
deltat = 0.02;
function affect!(integrator)
integrator.p[3] = rhseval(integrator.t)
end
cb = PeriodicCallback(affect!,deltat)
function rhseval(x)
return rhs[floor(Int,x/deltat)+1]
end
function sdof!(du,u,p,t)
omg = p[1]
zeta = p[2]
du[1] = u[2]
du[2] = - omg^2 * u[1] - 2zeta * omg * u[2] - p[3]
end
disp0 = 0;
velo0 = 0;
u0 = [disp0, velo0];
tspan = (0.0,0.1);
p = [2pi,0.02,0];
prob = ODEProblem(sdof!,u0,tspan,p,callback = cb)
sol = solve(prob,abstol = 1e-8, reltol = 1e-8,saveat=0.02)
but the obtained results are unsatisfactory.
Is there other way, not to interpolate g(t), but to use callback, i.e PeriodicCallback, DiscreteCallback?