CODEWITHSUNDEEP
pythonode
Faghetti
Faghetti

Reputation: 21

How do I make Python solve this Second-order nonlinear ODE?

I am currently working through some Lagrangian Mechanics, which is going pretty well for the most part, but I decided to try out a computer exercise and I can't figure out how to make Python solve these sorts of equations.

The Second-order nonlinear ODE I am trying to solve is:

enter image description here

with the intial value of y being y0 = 0.

Solving a normal, first order equation is easy as you just create a function set it equal to something and then use odeint. No problem. But now I got a second order equation and I do not know how to tell odeint or the function that what I'm trying to find is something differentiated twice.

I have tried to use various amounts of integration and differentiation functions (symbolic and otherwise) and said stuff like:

y' differentiated = RHS

y' = RHS integrated with respect to t

and so on and so forth, but nothing seems to work. Either it tells me that my function can't be called or that something doesn't have an attribute or that there is some invalid limits and so on and so forth.

Here is the code I have right now (barebone) where I haven't tried any sort of trickery to get it to work and it simply returns the result of the function if it was first order instead:

#Setting up the derivative
def model(y,t):
    M = g = R = 1
    m = 0.7

    k = g/(R*(M+m))
    return k*(m-M*np.sin(y))

#Initial condition
y0 = 0

#Time interval
t = np.linspace(0,20,1000)

#solve ODE
y = odeint(model,y0,t)

The expected result should be oscillating and crud, but of course this does not since it is wrong. I am admittedly pretty new to Python, and very bad at coding in general, so can anyone help me out here?

Upvotes: 2

Views: 801

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

You need to transform your second order equation into a first order system with 2 components

def model(u,t):
    y, v = u
    M = g = R = 1
    m = 0.7

    k = g/(R*(M+m))
    return [ v, k*(m-M*np.sin(y)) ]

Then your solution will also have two components per time index, so that you have to extract y as the first of them.

Upvotes: 2

Related Questions