Double pendulum problem

[Arati Bhattu]

Hello

I am trying to model a 2D double pendulum problem with the following geometry and initial conditions. However, the solution does not match the analytical results. I would like to verify whether my inputs are correct.

Also, if I plan to make one of the links flexible, what would be an efficient way to implement it?

The code is included below.

....................................................................................................................................................

import pychrono as chrono
import pychrono.irrlicht as chronoirr

# Create system

sys = chrono.ChSystemNSC()
sys.SetGravitationalAcceleration(chrono.ChVector3d(0, -9.81, 0))

# 1. Ground body (fixed)
ground = chrono.ChBody()
ground.SetFixed(True)
ground.SetPos(chrono.ChVector3d(0, 0, 0))
sys.Add(ground)

# 2. First pendulum body
pend1 = chrono.ChBody()
pend1.SetMass(0.79)
pend1.SetInertiaXX(chrono.ChVector3d(0.01, 0.01, 0.0658))
pend1.SetPos(chrono.ChVector3d(0, -0.5, 0))  # CG in center of bar
pend1.SetAngVelParent(chrono.ChVector3d(0, 0, 10))
sys.Add(pend1)

# Visual: 1m tall, half-dim = 0.5, shift visual down by 0.5
shape1 = chrono.ChVisualShapeBox(chrono.ChVector3d(0.01,1, 0.01))
shape1.SetColor(chrono.ChColor(0.5, 0.5, 0.8))
pend1.AddVisualShape(shape1)

# 3. Second pendulum body
pend2 = chrono.ChBody()
pend2.SetMass(0.3)
pend2.SetInertiaXX(chrono.ChVector3d(0.01, 0.01, 0.0250))
pend2.SetPos(chrono.ChVector3d(0, -1.5, 0))  # CG in center of second bar
pend2.SetAngVelParent(chrono.ChVector3d(0, 0, 5))
sys.Add(pend2)

# Visual for pend2
shape2 = chrono.ChVisualShapeBox(chrono.ChVector3d(0.01, 1, 0.01))
shape2.SetColor(chrono.ChColor(0.8, 0.5, 0.5))
pend2.AddVisualShape(shape2)

# 4. Revolute joint between ground and pend1 at origin
joint1 = chrono.ChLinkRevolute()
frame1 = chrono.ChFramed(chrono.ChVector3d(0, 0, 0))
joint1.Initialize(ground, pend1, frame1)
sys.Add(joint1)

# 5. Revolute joint between pend1 and pend2 at (0, -1, 0)
joint2 = chrono.ChLinkRevolute()
frame2 = chrono.ChFramed(chrono.ChVector3d(0, -1, 0))
joint2.Initialize(pend1, pend2, frame2)
sys.Add(joint2)

# 6. Visualization setup
vis = chronoirr.ChVisualSystemIrrlicht()
vis.AttachSystem(sys)
vis.SetWindowSize(1024, 768)
vis.SetWindowTitle("Double Pendulum")
vis.Initialize()
vis.AddSkyBox()
vis.AddCamera(chrono.ChVector3d(0.5, 1.5, 1.0))  # Camera looking in 3D
vis.AddTypicalLights()

# Simulation loop
while vis.Run():
    vis.BeginScene()
    vis.Render()
    vis.EndScene()
    sys.DoStepDynamics(0.005)

[Dario Mangoni]

Can you show the discrepancies that you are observing in the simulation?
Are you aware about the chaotic behaviour of the double pendulum? Double_pendulum#Chaotic_motion
A slight change in the initial conditions, little approximations, etc can lead to totally different results and this is somehow expected.

Dario

[Arati Bhattu]

Hi Dario,

Thank you for responding.

I have solved the same problem numerically using the geometry and initial conditions provided earlier. The expected output shows that the first link (pivoted to the ground) completes approximately one full revolution in a 0.8-second interval.

As I am still relatively new to Chrono, there is a possibility that I may not be providing the inputs as intended.

Best,

Arati

[Arati Bhattu]

Hi Dario,

Could you kindly look into the code? I also have tried running the simple simulation with single compound pendulum with initial velocity. I have figured out the 'pend1.SetAngVelParent(chrono.ChVector3d(0, 0, wz))' is not working the way I want to input it. For single pendulum problem with initial theta, the answer matches exactly with analytical solution. However with initial velocity (dtheta0) the frequency matches but amplitude. I am plotting X displacement of CG of pendulum . Here is the code for your reference and comparison plots

Mass =0.79 kg

length  =1 m

Initial velocity (dtheta0) =0.01 rad/sec

import pychrono as chrono
import pychrono.irrlicht as chronoirr

import numpy as np
from scipy.io import savemat
# ------------------------------------------------------------------------------
# 1. Create the Chrono physical system
# ------------------------------------------------------------------------------

sys = chrono.ChSystemSMC()
sys.SetGravitationalAcceleration(chrono.ChVector3d(0, -9.81, 0)) # Global gravity

# ------------------------------------------------------------------------------
# 2. Create the ground body (fixed, at global origin)
# ------------------------------------------------------------------------------

ground = chrono.ChBody()
ground.SetFixed(True)
ground.SetPos(chrono.ChVector3d(0, 0, 0))  # Explicit in global frame
sys.Add(ground)

# ------------------------------------------------------------------------------
# 3. Create the pendulum body (initial position, orientation, velocity all in global)
# ------------------------------------------------------------------------------

pend1 = chrono.ChBody()

# Set global position of pendulum CG
pend1.SetPos(chrono.ChVector3d(0, -0.5, 0))  # 0.5m below joint in Y (gravity) direction

# Set mass and inertia
pend1.SetMass(0.79)
pend1.SetInertiaXX(chrono.ChVector3d(0.01, 0.01, 0.0658))  # principal inertias
# Put dummy inertia in X and Y direction. Otherwise it gives Inf solution

# Set angular velocity in global frame (around Z-axis)
pend1.SetAngVelParent(chrono.ChVector3d(0, 0, 0.01))  # radians/sec in global Z
sys.Add(pend1)

# ------------------------------------------------------------------------------
# 4. Revolute joint between ground and pend1 (at global origin, around global Z)
# ------------------------------------------------------------------------------

# Frame at joint location (global coordinates)
joint_frame = chrono.ChFramed(chrono.ChVector3d(0, 0, 0))  # Revolute joint at origin
joint1 = chrono.ChLinkRevolute()
joint1.Initialize(ground, pend1, joint_frame)

sys.Add(joint1)

# Simulation parameters
end_time = 50.0  # seconds
step_size = 0.01
time = 0.0

# Initialize storage lists
time_data = []
pend1_x, pend1_y, pend1_z = [], [], []
pend2_x, pend2_y, pend2_z = [], [], []

# Run the simulation
while time < end_time:
    sys.DoStepDynamics(step_size)

    # Get positions
    pos1 = pend1.GetPos()

    # Append data to lists
    time_data.append(time)

    pend1_x.append(pos1.x)
    pend1_y.append(pos1.y)
    pend1_z.append(pos1.z)

    time += step_size

# Prepare data dictionary
mat_data = {
    'time': np.array(time_data, dtype=np.float64),
    'pend1_x': np.array(pend1_x, dtype=np.float64),
    'pend1_y': np.array(pend1_y, dtype=np.float64),
    'pend1_z': np.array(pend1_z, dtype=np.float64),
    'pend2_x': np.array(pend2_x, dtype=np.float64),
    'pend2_y': np.array(pend2_y, dtype=np.float64),
    'pend2_z': np.array(pend2_z, dtype=np.float64),
}

# Save to .mat file
savemat('trial', mat_data)

Comparison with analytical solution

[Dario Mangoni]

Can you try setting the initial linear velocity instead?

Do you notice any difference? 

-------- Messaggio originale --------

Da: Arati Bhattu <aa...@rice.edu>

Data: 05/06/25 01:36 (GMT+01:00)

A: ProjectChrono <projec...@googlegroups.com>

Oggetto: [chrono] Re: Double pendulum problem

--
You received this message because you are subscribed to the Google Groups "ProjectChrono" group.
To unsubscribe from this group and stop receiving emails from it, send an email to projectchron...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/projectchrono/264876b6-a1c3-44fa-b0c7-f8ef33f18aaan%40googlegroups.com.

[Dario Mangoni]

Hi Arati,
two ways to achieve proper results:

• setting initial velocities: in this case you need to properly set the initial velocities, that means both linear and angular; in fact the COG is not only rotating but also translating
      body1->SetPosDt({angvel_init * 0.5, 0, 0});
      body1->SetAngVelLocal({0,0,angvel_init});
• temporarily adding a speed-imposing motor and then removing it immediately after a DoAssembly

Here is my result and my CPP code (sorry I don't play much on PyChrono).

#include "chrono/physics/ChSystemNSC.h"
#include "chrono/physics/ChBodyEasy.h"
#include "chrono/physics/ChLinkMotorRotationSpeed.h"
#include "chrono/core/ChRealtimeStep.h"

#include "chrono_irrlicht/ChVisualSystemIrrlicht.h"

using namespace chrono;
using namespace chrono::irrlicht;

// -----------------------------------------------------------------------------

int main(int argc, char* argv[]) {
    std::cout << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << std::endl;

    // Create a Chrono physical system
    ChSystemNSC sys;

    double angvel_init = 0.01;

    // Create pendulum bodies
    auto body0 = chrono_types::make_shared<ChBody>();
    body0->SetFixed(true);
    sys.Add(body0);

    auto body1 = chrono_types::make_shared<ChBodyEasyBox>(0.1, 1, 0.1,  // x,y,z size
                                                            1,        // density
                                                            true,     // visualization?
                                                            false);   // collision?
    body1->SetPos(ChVector3d(0, -0.5, 0));
    body1->SetInertiaXX({0.01, 0.01, 0.0658});
    //body1->SetPosDt({angvel_init * 0.5, 0, 0});
    //body1->SetAngVelLocal({0,0,angvel_init});
    body1->SetMass(0.79);
    sys.Add(body1);

   

    // Create a spherical joint
    auto rev_joint = chrono_types::make_shared<ChLinkRevolute>();
    rev_joint->Initialize(body1, body0, ChFramed());
    sys.AddLink(rev_joint);

    auto mot = chrono_types::make_shared<ChLinkMotorRotationSpeed>();
    mot->SetMotorFunction(chrono_types::make_shared<ChFunctionConst>(angvel_init));
    mot->Initialize(body1, body0, ChFramed());
    sys.AddLink(mot);


    // Create a second spherical joint
    // Create the Irrlicht visualization system
    auto vis = chrono_types::make_shared<ChVisualSystemIrrlicht>();
    vis->AttachSystem(&sys);
    vis->SetWindowSize(800, 600);
    vis->SetWindowTitle("A simple pendulum example");
    vis->Initialize();
    vis->AddLogo();
    vis->AddSkyBox();
    vis->AddCamera(ChVector3d(0, 1, -1));
    vis->AddTypicalLights();


    // Simulation loop
    double timestep = 0.001;
    ChRealtimeStepTimer realtime_timer;

    sys.DoAssembly(AssemblyAnalysis::FULL);
    mot->SetDisabled(true);


    while (vis->Run()) {
        vis->BeginScene();

        vis->Render();
       
        tools::drawGrid(vis.get(), 2, 2, 20, 20, ChCoordsys<>(ChVector3d(0, -20, 0), QuatFromAngleX(CH_PI_2)),
                        ChColor(0.3f, 0.5f, 0.5f), true);


        vis->EndScene();
        std::cout << sys.GetChTime() << ", " << body1->GetPos().x() << "; ";



        sys.DoStepDynamics(timestep);


        realtime_timer.Spin(timestep);
    }

    return 0;
}

[Arati Bhattu]

Hi Dario,

Thank you for your reply. The solution works well and matches closely in Python using the two approaches you described. I will definitely keep in mind the importance of providing consistent boundary conditions.

I have just one more request: could you please check your solution for a large initial angle or a large initial angular velocity? (For example, I used an initial angular velocity of 5 rad/sec.) I tried running the same code with the modified nonlinear analytical equation—removing the small-angle assumption. The solutions match well for the first 20 seconds, but after that, differences start to appear.

Thank you again for your help!

Best,

Arati

You received this message because you are subscribed to a topic in the Google Groups "ProjectChrono" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/projectchrono/eZBeFTKyJjY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to projectchron...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/projectchrono/ca2240a7-f99a-4063-84e4-df7fcced17adn%40googlegroups.com.

[Arati Bhattu]

Hi Dario,

I just wanted to share that I realized the solution tolerance used when solving the ODE makes a significant difference. The discrepancies I observed earlier were due to differences in the tolerance values between the implementations.

Thanks,

Arati