# Pendulum using Crank-Nicolson

##### 3 Comments

### Answers (1)

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.

##### 4 Comments

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');

### See Also

### Products

### Community Treasure Hunt

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

Start Hunting!