Hi, new to Julia; the instructions on how to use “SecondOrderODE” on the julia Diff Eq package website are very incomplete. What is wrong with the very simple program below? Thanks.
BTW, I’m just trying to model the Simple Harmonic Oscillator: u" = -9 u .
using DifferentialEquations
function OneDimSpring!(dv,v,u,p,t)
du = v
dv = -9u
end
x0 = [2.0]
v0 = [0.0]
tspan = (0.0,30.0)
prob = SecondOrderODEProblem(OneDimSpring!,v0,x0,tspan)
sol = solve(prob)
I get the error message on solve: "MethodError: no method matching similar (::Float64, :: {TypeFloat64}). Putting specific Diff Eq solution commands in the “prob” statement doesn’t fix it. (As an alternative to du = v, I also tried dv = u. That didn’t work either.)
I’m trying the in-place solution method; I tried function OneDimSpring!(ddu,du,u,p,t), but I couldn’t figure out the right syntaxx for that either.
In your function, du isn’t (and doesn’t need to be) defined. The assumption for SecondOrderODEProblems is that du = v, so it isn’t supposed to be written in your dynamics function.
Reformatting for readability and fixing some minor syntax, here is your original function.
function OneDimSpring!(du,u,p,t)
du[1] = u[2]
du[2] = -9u[1]
end
This is already in the form of a 1st order ODE and can be solved with ODEProblem
function spring_manually_reduced()
x0 = [2.0, 0.0]
tspan = (0.0,30.0)
prob = ODEProblem(OneDimSpring!,x0,tspan)
sol = solve(prob)
end
SecondOrderODEProblem does this reduction for you so you only need to specify the function to be double integrated.
function OneDimSpring2!(ddu, du, u, p, t)
ddu[1] = -9u[1]
end
function spring2()
du0 = [0.0]
u0 = [2.0]
tspan = (0.0,30.0)
prob = SecondOrderODEProblem(OneDimSpring2!,du0,u0,tspan)
sol = solve(prob)
end
using DifferentialEquations
function OneDimSpring!(dv,v,u,p,t)
dv .= -9u
end
x0 = [2.0]
v0 = [0.0]
tspan = (0.0,30.0)
prob = SecondOrderODEProblem(OneDimSpring!,v0,x0,tspan)
sol = solve(prob)
Thanks for your help, I got the thing to compile & plot now, but the solution is wrong – a straight line instead of a sinusoid. This is how I streamlined it (still to get u" = -9u) :
What’s up? (Thanks.)
function OneDimSpring(dv,v,u,p,t)
dv = -9u
end
x0 = [2.0]
v0 = [1.0]
tspan = (0.0,30.0)
prob = SecondOrderODEProblem(OneDimSpring,v0,x0,tspan)
sol = solve(prob)
.= is because the variables are vectors, so you need to apply to each element of the vector (which is only one, but still, vectors are vectors and scalars are scalars, unlike, say in MATLAB.) The dot syntax in general in julia means “broadcast this to every element in the collection”. So A .= B means “set each element of A to its corresponding element in B“. It’s something you’ll see a lot and it’s super useful.
The other option is to do what @contradict did and index into your variables directly.
Here is the relevant section of the julia docs on the dot syntax. Even though the manual is pretty long, it’s very clear and well written IMO and definitely worth reading.
Thanks guys, that makes sense and it all works now, I’m happy. I will surely have more questions.
For now, I’m wondering about the simplest way to integrate a Diff Eq backwards in time, so that my “initial conditions” are applied to the highest time value (or perhaps somewhere in the middle of the time range).
Would something like tspan = (0.0, -10.0) be enough to integrate from 0 BACKWARDS to -10?
And how about tspan = (10.0, 0.1), to go from 10.0 back to 0.1?
Or better yet, a way of applying my initial conditions to, say, t = 1.0, while integrating backwards and forward from that point, say, tspan = (0.0,10.0) ?
I’m sure there’s a way of specifying backwards timesteps, like perhaps dt = -0.1, but wouldn’t that prevent the solver from using variable timesteps in the case of a difficult equation?
I’m asking this because my physics problem blows up at t=0, and it doesn’t make sense to define initial conditions there. And negative time has no meaning here.
I’ve never actually had a reason to try integrating backwards in time, but I think it might be something that’s handled. I’m not at a computer right now, but you can give it a shot with (10.0, 0.0). If it fails, then I guess that isn’t the way to do it.
Specifying the initial conditions anywhere near t=0 don’t work, because it’s Bessel functions (in the real equations), which either blow up or go to zero near there. (And t=0 is the Big Bang, so I can’t start the calculation at negative t, which is unphysical.)
I really need to see the behavior close to t → 0, without starting the problem there.
Setting tspan = (10.0,0.0) actually did make the intitial conditions be applied at t =10 instead of t=0, thanks. (The plot is also backwards, a minor inconvenience, but I should be able to fix that.)
BTW, some of my 2nd-order equations have non-homogeneous terms (nonzero right-hand-sides). But they are simpler as 4th-order equations. Does anyone know if Julia has a “FourthOrderOde” command, as an extension of its “SecondOrderODE” command?
Generally there aren’t any algorithms out there specialized for solving 4th order ODEs, so I’d just treat it through order lowering. You can follow the README example of ModelingToolkit.jl as a nice quick way to do lowering: