Plot Lagrangian solutions

Question

I know R an Python, could you help me find a simple technique to plot the solutions of a system of differential equations?(In this case a Lagrangian)

For example, I have this system:

$$ kx=mx''+ml\theta''\cos(\theta)-ml\theta'^2\sin(\theta)\\ -mlx'\theta'\sin(\theta)-mgl\sin(\theta)+k\frac{l}{2}(l+\frac{l}{2}\cos(\theta))\sin(\theta)=ml^2\theta''+mlx''\cos(\theta)-mlx'\theta'\sin(\theta)+\frac{m}{8}l^2\theta'' $$ Where x and $\theta$ are my variables, and the others letters are parameters.

This is the solution I found from an exercise I'm doing, but I don't know if I'm correct because I don't have the solutions, so, maybe visualising the graph I understand.

Thank you.

Answer 1

If you can figure out Scipy and matplotlib, then you can solve and plot ODEs with the least amount of effort. Use the integrate function from Scipy to solve the system and matplotlib to plot it. However, it's also possible to achieve what you want without installing anything new, so I'll describe that procedure here. Regardless of which method you use, we still have to put the DEs into a more suitable form. Specifically, we must convert the system of two coupled second-order equations into four coupled first-order equations (you can always do this.)

The first step is to solve each of the equations for $x''$ and $\theta''$, note that we can isolate $\theta ''$ by solving for $x''$ in the first equation and eliminating it from the second. This will give $\theta ''$ and $x''$ as (messy) functions of $x, \theta, x'$ and $\theta'$. Now we make the substitutions $y_{1} = x, y_{2} = \theta, y_{3} = x', y_{4} = \theta '$. Note that we have $y_{1}' = y_{3}$ and $y_{2}' = y_{4}$. Furthermore, $y_{3}' = x''$ and $y_{4}' = \theta ''$. But since we just solved for $x''$ and $\theta''$, we can use those here to get $y_{3}'$ and $y_{4}'$ as functions of $y_{1-4}$. The net result is that we have a four dimensional vector $\vec{y}$ whose first derivative $\vec{y}\,'$ we know. See Example $1$ here for an example of this procedure being carried out explicitly.

Now that we have a first-order system, we can use an integration scheme to solve it. Fourth-order Runge-Kutta is easy to use; here is a quick Python script I wrote that implements it:

def mult(number, vector): # Multiplies "number" by each of the elements of "vector"
    vec = []
    for x in range(0, len(vector)):
        vec.append(number*vector[x])
    return vec


def add(vectors): # Adds list of vectors to eachother component-wise
    vec = []
    for x in range(0, len(vectors[1])):
        tot = 0.0
        for y in range(0, len(vectors)):
            tot = tot + vectors[y][x]
        vec.append(tot)
    return vec


def make_string(t, y): # Formats scalar t and vector y to be suitable for writing to file 
    s = ''
    s = s + str(t) + ','

    for x in range(0, len(y)):
        if x < len(y) - 1:
            s = s + str(y[x]) + ','
        else:
            s = s + str(y[x]) +'\n'
    return s


def derivatives(t, y): # Returns the derivatives of y evaluated at a particular time and position
    y1_prime = y[2]
    y2_prime = y[3] 
    y3_prime = -y[0] # x''
    y4_prime = -y[1] # theta''
    return [y1_prime, y2_prime, y3_prime, y4_prime]


def RK4(): # Core Runge-Kutta procedure
    t = 0.0 # Initial time
    y = [0.0, 0.0, 0.1, 0.1] # Initial positions 
    dt = 0.01 # Time_step
    t_max = 10.0 # Time to integrate

    f = open('file.txt', 'w') # Output file
    f.write('t,x,theta,x_prime,theta_prime\n') # Column names
    f.write(make_string(t, y)) # Write initial values

    while t < t_max:
        k1 = derivatives(t, y)
        k2 = derivatives(t + 0.5*dt, add([y,mult(0.5*dt,k1)]))
        k3 = derivatives(t + 0.5*dt, add([y,mult(0.5*dt,k2)]))
        k4 = derivatives(t + dt, add([y,mult(dt,k3)]))

        increment = mult(1.0/6.0, add([k1, mult(2.0,k2), mult(2.0,k3), k4]))

        y = add([y, mult(dt,increment)])
        t = t + dt

        f.write(make_string(t, y))
    f.close()
    return

This example solve two simultaneous simple harmonic oscillator problem simultaneously. For your plot, you'll need to modify the expressions for the derivatives in the function derivatives to what I described earlier. Furthermore, you'll have to set the initial conditions. Running RK4() generates a file called file.txt which contains a list of points with column names. Next, enter the following commands into $R$ to make a plot (you'll have to set the working directory to point to the location of file.txt):

data <- read.csv(file="file.txt",head=TRUE,sep=",")

plot(data$t, data$x) or plot(data$t, data$theta)

This should generate the required plots.