Simple 2nd Order Differential Equation syntax?

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)

Also, jonniedie, what’s with the “dot-equals” ?

.=

??

You need to either index or broadcast into dv. Check out the examples above to see the two options. What you’re doing here is overwriting dv

.= 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.

Oh, by the way, welcome!

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.

Thanks again…

1 Like

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.

If it’s t=0.0 that’s giving you the problem, can’t you start it at 0.1? so tspan = (0.1, 10.0)

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:

1 Like

OK, thanks, I’ll look into it!