How has the author mathematically modeled a double pendulum in this study?

Question

I'm looking for project ideas for my final year dissertation for a computing degree. Here is the study (A double pendulum model of tennis strokes. Rod Cross. Uni of Sydney, 2006) - http://www.physics.usyd.edu.au/~cross/PUBLICATIONS/49.%20TennisDPend.pdf. My matlab files and spreadsheet files are linked at the bottom of the post.

My goal was to investigate a link between arm length/mass and the angular velocity/swingspeed of a forehand, mapped by a pendulum .

I've looked into this before, but I wasn't able to solve the calculus. I'm a fairly strong coder but I'm terrible with calculus. I've coded the pendulum in an executable matlab code, but the angular velocity peaks too early in time, making the pendulum swing too fast, suggesting there is an issue with the formulas. Is there any other pendulums that I could use, that are solvable.

I've put all ODE's into a spreadsheet. The author doest state which numerical method was used. So I implemented the formulas using the RK4 method, to try and investigate the formulas, and why the velocity is rising too early. It's not made clear in the paper how L2, M2, h2 and I2 are calculated. It's also not mentioned whether the hand was included in the calculations for the lower pendulum, which would really complicate it, as h2 and the moment of inertia could differ, depending on the position of the hand on the grip, the mass distribution of the hand and the size of the hand.

I also tried to replicate the results from this study , which is based on the same equations from Rod Cross' study to see if that helped - https://arxiv.org/pdf/1505.01916.pdf . I had the exact same problem of the chart looking similar but rising too high, too early. .https://arxiv.org/pdf/1505.01916.pdf. I emailed the authors of both studies to check what method/timestep they used, but never heard back. I also asked whether the hand was included.

Matlab files are here - https://ufile.io/f/s3jmi

Spreadsheet is here. - https://drive.google.com/file/d/1uzjtXV2bweHj7ADU47lX3NaNEX7mZPEK/view?usp=sharing

Answer 1

Seems like the author in the first paper is using standard Newtonian analysis.

(EDIT: If you get sick of pendula, too, you can always look at masses on springs.)

Not sure if I'll get downvoted for this move, but I will paste some code below that solves and plots (in real time) the double pendulum using the RK4 method. Problem is, it's QB64 code, which you're welcome to attain and run yourself - but I'm sure you want the equations (just for comparison to your own) more than the program. Here it all is in context:

DEFSNG A-Z
pi = 3.141593
r = 100
a1 = -pi / 2
a2 = -pi / 2
h = 0.001
m = 1000
g = 1000
sw = 640
sh = 480

DIM p1 AS LONG
p1 = _NEWIMAGE(sw, sh, 12)
SCREEN _NEWIMAGE(sw, sh, 12)

DO
    a1p = (6 / (m * r ^ 2)) * (2 * m * v1 - 3 * COS(a1 - a2) * m * v2) / (16 - 9 * (COS(a1 - a2)) ^ 2)
    a1k1 = a1p
    a1k2 = a1p + h * a1k1 / 2
    a1k3 = a1p + h * a1k2 / 2
    a1k4 = a1p + h * a1k3
    a2p = (6 / (m * r ^ 2)) * (8 * m * v2 - 3 * COS(a1 - a2) * m * v1) / (16 - 9 * (COS(a1 - a2)) ^ 2)
    a2k1 = a2p
    a2k2 = a2p + h * a2k1 / 2
    a2k3 = a2p + h * a2k2 / 2
    a2k4 = a2p + h * a2k3
    na1 = a1 + h * (a1k1 + 2 * a1k2 + 2 * a1k3 + a1k4) / 6
    na2 = a2 + h * (a2k1 + 2 * a2k2 + 2 * a2k3 + a2k4) / 6
    v1p = -0.5 * r ^ 2 * (a1p * a2p * SIN(a1 - a2) + 3 * g * SIN(a1) / r)
    v1k1 = v1p
    v1k2 = v1p + h * v1k1 / 2
    v1k3 = v1p + h * v1k2 / 2
    v1k4 = v1p + h * v1k3
    v2p = -0.5 * r ^ 2 * (-a1p * a2p * SIN(a1 - a2) + g * SIN(a2) / r)
    v2k1 = v2p
    v2k2 = v2p + h * v2k1 / 2
    v2k3 = v2p + h * v2k2 / 2
    v2k4 = v2p + h * v2k3
    nv1 = v1 + h * (v1k1 + 2 * v1k2 + 2 * v1k3 + v1k4) / 6
    nv2 = v2 + h * (v2k1 + 2 * v2k2 + 2 * v2k3 + v2k4) / 6
    a1 = na1
    a2 = na2
    v1 = nv1
    v2 = nv2

    _DEST p1
    PSET (sw / 2 + r * COS(a1 - pi / 2) + r * COS(a2 - pi / 2), sh / 2 - r * SIN(a1 - pi / 2) - r * SIN(a2 - pi / 2)), 12
    _DEST 0
    _PUTIMAGE , p1

    LINE (sw / 2, sh / 2)-STEP(r * COS(a1 - pi / 2), -r * SIN(a1 - pi / 2))
    CIRCLE STEP(0, 0), 5
    LINE -STEP(r * COS(a2 - pi / 2), -r * SIN(a2 - pi / 2))
    CIRCLE STEP(0, 0), 5

    _LIMIT 200
    _DISPLAY
LOOP UNTIL _KEYHIT = 27
SYSTEM