# Pendulum using Crank-Nicolson

7 views (last 30 days)
FW on 6 May 2019
Commented: Umar on 12 Jul 2024
Hi everybody,
I'm relatively new to MatLab, and to numerical analysis in general. I have this problem I have to solve, using the Crank-Nicolson method, and I don't have the slightest idea how to start my code...
I'm trying to solve two non-linear forced pendulum equation, which I adimensionalized as follows :
Can anyone give me a hint on how to start my process?
Thanks!
Kena
##### 3 CommentsShow 1 older commentHide 1 older comment
David Goodmanson on 12 Jul 2024
Edited: David Goodmanson on 12 Jul 2024
Is there a reference available for this set of equations?
Umar on 12 Jul 2024
Hi David,
The reference was already provided by Shresth, please correct me if I am wrong.

Umar on 12 Jul 2024

Hi FW,

You asked, Can anyone give me a hint on how to start my process?

Well, I can walk you through it to help you understand and start your journey which will help you implement your code. First I defined various constants related to the physical properties of the pendulum system, such as acceleration due to gravity, length of the pendulum, mass, amplitude of the forcing function, natural frequency, and damping coefficient listed below.

Constants Initialization

• g = 9.81: Represents the acceleration due to gravity.
• L = 1.0: Denotes the length of the pendulum.
• m = 1.0: Indicates the mass of the pendulum.
• A = 0.5: Represents the amplitude of the forcing function.
• omega_0 = 2.0: Signifies the natural frequency of the pendulum.
• beta = 0.1: Represents the damping coefficient.

Time Parameters

• dt = 0.01: Defines the time step for the simulation.
• T = 10: Specifies the total simulation time.
• t = 0:dt:T: Creates a time vector from 0 to T with increments of dt

Initial Conditions

• theta and omega are initialized as arrays of zeros with the same size as the time vector t.
• theta(1) = 0.1: Sets the initial angle of the pendulum.
• omega(1) = 0: Sets the initial angular velocity of the pendulum.

Then, I applied Crank-Nicolson method which is used to numerically solve the differential equations governing the motion of the damped pendulum system. It involves updating the values of theta and omega at each time step based on the previous values and the system dynamics. I implemented the following equations are used in the loop:  theta(i) = (theta(i-1) + dt*omega(i-1)) / (1 + 0.5*dt^2*(omega(i-1)^2 - g/L*cos(theta(i-1)) + A*cos(omega_0*t(i-1))));   omega(i) = (omega(i-1) - 0.5*dt*(g/L*sin(theta(i-1)) + beta*omega(i-1) - A*sin(omega_0*t(i-1)))) / (1 + 0.5*dt*beta);   Afterwards, I draw conclusion by plotting the simulation results in two subplots: 1. Subplot 1: * Plots theta against time. * Sets the x-axis label to 'Time' and the y-axis label to 'Theta'. * Titles the plot as 'Theta vs Time'. 2. Subplot 2: * Plots omega against time. * Sets the x-axis label to 'Time' and the y-axis label to 'Omega'. * Titles the plot as 'Omega vs Time'.

Please see detailed code below as reference along with attached plots.

>> % Constants g = 9.81; % acceleration due to gravity L = 1.0; % length of the pendulum m = 1.0; % mass of the pendulum A = 0.5; % amplitude of the forcing function omega_0 = 2.0; % natural frequency of the pendulum beta = 0.1; % damping coefficient

% Time parameters dt = 0.01; % time step T = 10; % total simulation time t = 0:dt:T; % time vector

% Initial conditions theta = zeros(size(t)); omega = zeros(size(t)); theta(1) = 0.1; % initial angle omega(1) = 0; % initial angular velocity

% Crank-Nicolson method for i = 2:length(t) theta(i) = (theta(i-1) + dt*omega(i-1)) / (1 + 0.5*dt^2*(omega(i-1)^2 - g/L*cos(theta(i-1)) + A*cos(omega_0*t(i-1)))); omega(i) = (omega(i-1) - 0.5*dt*(g/L*sin(theta(i-1)) + beta*omega(i-1) - A*sin(omega_0*t(i-1)))) / (1 + 0.5*dt*beta); end

% Plotting results figure; subplot(2,1,1); plot(t, theta, 'b'); xlabel('Time'); ylabel('Theta'); title('Theta vs Time');

subplot(2,1,2); plot(t, omega, 'r'); xlabel('Time'); ylabel('Omega'); title('Omega vs Time');

Attached code and plots

By understanding each step of the code, you can grasp how the Crank-Nicolson method is applied to simulate the behavior of a damped pendulum system over time.

Umar on 12 Jul 2024

Please ignore my most recent comment. However, I made some tweaks to the code, please see below and provide your suggestions if it looks okay to you. When you are implementing complex equations in code, it can create bugs. So, it is always good idea to have second set of eyes.

% Constants

g = 9.81; % acceleration due to gravity

L = 1.0; % length of the pendulum

m = 1.0; % mass of the pendulum

A = 0.5; % amplitude of the forcing function

omega_0 = 2.0; % natural frequency of the pendulum

beta = 0.1; % damping coefficient

% Time parameters

dt = 0.01; % time step

T = 1000; % total simulation time

t = 0:dt:T; % time vector

% Initial conditions

theta = zeros(size(t));

omega = zeros(size(t));

theta(1) = 0.1; % initial angle

omega(1) = 0; % initial angular velocity

% Crank-Nicolson method

for i = 2:length(t)

`    theta(i) = (theta(i-1) + dt*omega(i-1)) / (1 + 0.5*dt^2*(omega(i-1)^2 - g/L*cos(theta(i-1)) + A*cos(omega_0*t(i-1))));`
```    omega(i) = (omega(i-1) - 0.5*dt*(g/L*sin(theta(i-1)) + beta*omega(i-1) - A*sin(omega_0*t(i-1))) / (1 + 0.5*dt*beta));
end```

% Plotting results

figure;

subplot(2,1,1);

plot(t, theta, 'b');

xlabel('Time');

ylabel('Theta');

title('Theta vs Time');

subplot(2,1,2);

plot(t, omega, 'r');

xlabel('Time');

ylabel('Omega');

title('Omega vs Time');

David Goodmanson on 12 Jul 2024
Edited: David Goodmanson on 12 Jul 2024
Yes the result is different from before, but I realized I don't know enough about the initial set of equations to come to any sensible conclusion and I modified a previous comment accordingly.

### Categories

Find more on Programming in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!