Differential Equation of a Pendulum

[Tygra]

Homework Statement

How to compute the angle of displacement

Relevant Equations

In Post

Hi everyone,

I am going through the book called "Numerical Methods for Engineers" and I am trying to solve for the angle of displacement of a pendulum. I am using MATLAB to do this. However, I am getting some strange results.

In the book the differential equation is:

Because this is a second order differential equation, I have converted it to two first order equations so that I can solve it using Euler's method.

So I get:

Where:

g = gravitational constant
l = Pendulum length

g = 9.81;
l = 0.2;

t = 0:0.1:50;
dt = t(2) - t(1);

z(1) = 0
theta(1) = 90;

for i = 1:length(t)-1
    theta(i+1) = z(i)*dt + theta(i)
    z(i+1) = -(g/l)*sind(theta(i))*dt + z(i)
  
end

This is the graph

This does not look correct to me. The amplitude of the graph should decrease as the time goes on correct?

Can someone help please?

 

[haruspex]

Any attempt to model physical processes on the basis of forces and accelerations suffers from cumulative errors in the integration. A steady increase in the system's energy is not unusual.
Reducing the time step will improve it but not cure it.

The amplitude of the graph should decrease as the time goes on correct?

No, because you have not modelled any losses. The true result would be constant amplitude.

The eventual runaway in your chart is where the angle reaches 180°, so the pendulum goes over the top, spins right round, and goes faster and faster.

A crude fix is to assess the total energy at each step and tweak the velocity or angle to keep it constant. That's not really satisfactory because you could not use that if modelling friction, for example.

To understand the details, consider the downward trajectory. You calculate the acceleration for the next interval, but the acceleration is reducing, so the acceleration you use is greater than the average acceleration will be.

In your code, you calculate the new position then calculate the new velocity (z) based on the acceleration at the old position.

Edit: See post #8
You would do better to base it on some weighted average of the old and new positions. Could be a simple average, or there might be a better weighting.
Whether it is better to average the positions then take the sines or take the average of the sines I'm not sure without more detailed analysis.

 

[vela]

One problem is that angles are measured in radians in math. In your numerical solution, you're using degrees. You'll need to put in a conversion factor, or more simply, use radians in the numerical solution.

Ideally, the amplitude of the oscillation shouldn't decrease because you haven't included any friction.

 

[haruspex]

One problem is that angles are measured in radians in math. In your numerical solution, you're using degrees. You'll need to put in a conversion factor, or more simply, use radians in the numerical solution.

I assume "sind" means the argument is in degrees. Other than the error accumulation the results look right.

Ideally, the amplitude of the oscillation shouldn't decrease because you haven't included any friction.

But it shouldn’t increase either.

 

[vela]

I assume "sind" means the argument is in degrees. Other than the error accumulation the results look right.

But $g/l$ in the differential equation isn't in (degrees/s)^2. Given the values the OP used, the period should be about 0.9 s.

But it shouldn’t increase either.

I was addressing the OP's belief the amplitude should decrease with time. I figured you already talked about why it increased in the numerical solution.

 

[haruspex]

But $g/l$ in the differential equation isn't in (degrees/s)^2. Given the values the OP used, the period should be about 0.9 s.

True, so that means the effective time units being used are not seconds.

 

[Tygra]

Thank you guys for your response. Considering what has been said how would I amend the equation to include friction/drag, so that the model behaves realistically?

 

[haruspex]

Thank you guys for your response. Considering what has been said how would I amend the equation to include friction/drag, so that the model behaves realistically?

I've thought of a better (and more natural) way. Just use more terms of the Taylor series.
$\omega(t+\Delta t)\approx\omega(t)+\alpha(t)\Delta t+\frac 12\dot\alpha(t)$. Similarly for $\theta$.
That generalises to the lossy case provided the nonconservative force is differentiable.

For the lossless case you can do even better by using the above for computing the new angle but then using energy conservation for the new angular velocity.

 

[kuruman]

I've thought of a better (and more natural) way. Just use more terms of the Taylor series.
$\omega(t+\Delta t)\approx\omega(t)+\alpha(t)\Delta t+\frac 12\dot\alpha(t)$. Similarly for $\theta$.

Should it not be
$\omega(t+\Delta t)\approx\omega(t)+\alpha(t)\Delta t+\frac 12\dot\alpha(t)(\Delta t)^2~$?

 

[haruspex]

Should it not be
$\omega(t+\Delta t)\approx\omega(t)+\alpha(t)\Delta t+\frac 12\dot\alpha(t)(\Delta t)^2~$?

Thanks for spotting it! Didn’t make it from brain to fingers.