Solve pendulum differential equation with for loop and no small angle approximation

Hi I am supposed to write a code that would solve a pendulum differential equation without the use of ode.

I cannot use the ode function in MATLAB and am supposed to use a for loop. I also cannot use the small angle approximation in the equation. This is where I have the issue because I am not sure how to incorporate into the code the sin function and have it run correctly. I someone could help me figure out where I have incorrect code that would be great.

This is the code that I have:

There error comes up in line 13 within the for loop

m=10;
l=1; %rope length
g=9.8; %g constant
syms x(t);
Dx= diff(x,t);
ode = l*diff(Dx)+g*sin(x(t))==0;
p=deg2rad(179);
cond=x(0)==p;
cond1=Dx(0)==0;
xSol(t)=dsolve(ode,[cond cond1]);
n=linspace(0, 10, 100);
for i=1:length(n)
    sol(i)=double(xSol(n(i)));
end
plot(n,sol);

[Bjorn Gustavsson]

From your description it seems as you are expected to implement an explicit Euler scheme for your equation...

Answer 1

[Milan Bansal]

Hi Alyson,

I understand that you are facing an error when trying to solve a pendulum differential equation without using MATLAB's ODE functions and incorporating the sine function directly. The issue arises from trying to assign values to "sol(i)" based on the evaluation of "xSol" at each point in n.

Solving a second-order differential equation representing a pendulum's motion without the small angle approximation makes the problem nonlinear due to the "sin(x(t))" term. The direct use of symbolic solutions "dsolve" is not allowed as per the constraints and a numerical solution with a for loop needs to be implemented.

To numerically solve the pendulum equation, discretize the problem using a simple numerical method like Euler's method.

1. The pendulum's equation of motion, not using the small angle approximation, is , where ( θ ) is the angle of the pendulum from the vertical. Convert this second-order ODE into two first-order ODEs.
2. Use a for loop to iterate over time steps and update the angle and angular velocity at each step.

% Pendulum parameters
m = 10; % mass (not used in the equations, but might be relevant for a more complex model)
l = 1; % rope length 
g = 9.8; % gravity

% Initial conditions
theta_0 = deg2rad(179); % initial angle in radians
omega_0 = 0; % initial angular velocity 

% Time parameters
T = 10; % total time 
dt = 0.1; % time step
time = 0:dt:T; % time vector

% Initialization
theta = zeros(size(time));   omega = zeros(size(time));
theta(1) = theta_0;  omega(1) = omega_0;

% Euler's method to solve the ODE
for i = 1:length(time)-1
    omega(i+1) = omega(i) - (g/l)*sin(theta(i))*dt;
    theta(i+1) = theta(i) + omega(i+1)*dt;
end

% Plot
plot(time, theta);
xlabel('Time (s)');
ylabel('Angle (rad)');

Hope this helps!