Pendulum using Crank-Nicolson

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

[Shreshth]

Hey,

you can refer to the below discussion to know more about the Crank-Nicolson applications using MATLAB.

http://www.mathworks.com/matlabcentral/answers/1653300-solve-the-partial-differential-equation-using-crank-nicolson-method

Thanks.

[David Goodmanson]

Is there a reference available for this set of equations?

[Umar]

Hi David,

The reference was already provided by Shresth, please correct me if I am wrong.

Answer 1

[Umar]

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.

[David Goodmanson]

Since the code was posted as an image rather than text that can be copied, I assumed the code in text above the image, which can be copied, was the same. I ran the code out to t = 100 sec which produces the following result.

[Umar]

Hi David,

Could you share the screenshot of my script along with the plot, I will appreciate that.

[Umar]

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]

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.