Differential Eqns in Mathematica help please

• Apr 8th 2010, 08:21 PM
scottie
Differential Eqns in Mathematica help please
Hey all.

I've got an assignment for a theoretical physics course which basically requires mathematica, and I've never been taught how to use it properly so some help would be good.

The task is graph out the motion of a double pendulum (one pendulum is attached to a point, the other attached at some arbitrary distance along the first pendulum)

So taking theta1 as the angle by which the first pendulum differs from pointing straight down, and theta2 as the angle by which the second pendulum differs from pointing straight down, I have got some equations of motion and plugged them into mathmatica in what I had hoped was the correct way, however for some initial conditions I get told that I am dividing by 0, or occasionally I get the comment "step size is effectively zero, singular or stiff system suspected"

Anyway, code is below, if anyone knows what I'm doing wrong, then that'd be awesome

Clear[t]
s = NDSolve[{theta1''[t] + Cos[theta1[t] - theta2[t]] theta2''[t] +
Sin[theta1[t] - theta2[t]] theta1'[t] theta2'[t] +
Sin[theta1[t]] == 0,
theta2''[t] +
Cos[theta1[t] - theta2[t]] theta1''[
t] - (theta1'[t]) (theta2'[t]) Sin[theta1[t] - theta2[t]] +
Sin[theta2[t]] == 0, theta1[0] == 2.0, theta1'[0] == 2.1,
theta2[0] == 1.0, theta2'[0] == 1.0001}, {theta1, theta2}, {t, 0,
10}]

Plot[Evaluate[{theta1[t] /. s, theta2[t] /. s}], {t, 0, 10},
PlotRange -> All]

Note - these initial conditions work, but for none working ones, an example is the case where both pendulums are pointing straight down, with no initial velocity, ie

Clear[t]
s = NDSolve[{theta1''[t] + Cos[theta1[t] - theta2[t]] theta2''[t] +
Sin[theta1[t] - theta2[t]] theta1'[t] theta2'[t] +
Sin[theta1[t]] == 0,
theta2''[t] +
Cos[theta1[t] - theta2[t]] theta1''[
t] - (theta1'[t]) (theta2'[t]) Sin[theta1[t] - theta2[t]] +
Sin[theta2[t]] == 0, theta1[0] == 0, theta1'[0] == 0,
theta2[0] == 0, theta2'[0] == 0}, {theta1, theta2}, {t, 0, 10}]

Plot[Evaluate[{theta1[t] /. s, theta2[t] /. s}], {t, 0, 10},
PlotRange -> All]

Thanks!

- Scott