Question

Solve the frictionless pendulum equation use Runge- kuta 4 or forward Euler

ly'' = -g*sin(y)

where g is gravitational acceleration and l is the length of the pendulum. The function y(t) represents the angle of the pendulum with respect to the vertial and y'(t0 the angular velocity. You will need to write the second-order equation as a system of two first-order equations, and you will need to write a function file that will evaluate this system of equations. Plot the solutions y1(t) and y2(t) on the same set of axes. Give an interpretation of y2(t). Set the length to be 1m, and use as initial conditions an initial displacement of pi/2 radians and a 0 rad/sec angular velocity

Answer 1

The Matlab scripts rk4sys.m, dydtSys.m, rk4SysSolver.m are posted below along with the screen shot of the generated plot. Run the script rk4SysSolver.m to generate the plot.

Also, note that the system of first order differential equations are solved using Runge Kutta 4 numerical technique.

rk4sys.m

function [t, yp] = rk4sys(dydt, tspan, y0, h)
% rk4sys: fourth-order Runge-Kutta for a system of ODEs
% [t, y] = rk4sys(dydt, tspan, y0, h): integrates
% a system of ODEs with fourth-order RK method

% input:
% dydt = name of the M-file that evaluates the ODEs
% tspan = [ti, tf]; initial and final times with output
% generated at interval of h, or
% = [t0 t1 ... tf]; specific times where solution output
% y0 = initial values of dependent variables
% h = step size

% output:
% t = vector of independent variable
% yp = vector of solution for dependent variables

ti = tspan(1);
tf = tspan(2);

t = ti: h: tf;
y(1, :) = y0;
yp(1, :) = y(1, :);

for i = 1: length(t)-1
    k1 = dydt( t(i), y(i,:) );
    ymid = y(i, :) + k1 * h/2;
    k2 = dydt( t(i) + h/2, ymid );
    ymid = y(i, :) + k2 * h/2;
    k3 = dydt( t(i) + h/2, ymid );
    yend = y(i, :) + k3 * h;
    k4 = dydt( t(i) + h, yend );
    phi = (k1 + 2 * (k2 + k3) + k4)/6;
    y(i + 1, :) = y(i, :) + phi * h;
    yp(i+1, :) = y(i+1, :);  
end
end

dydtSys.m

function rhs = dydtSys(t, y)
% governing ODE dydt = rhs
% dydt = dydt(t, y)

% input:
% x = independent variable
% y = dependent variable

% output:
% dydt = rhs of the ODE

g = 9.81;
l = 1   ;

dydt1 = y(2);
dydt2 = -g/l * sin( y(1) );
rhs = [dydt1, dydt2];

end

rk4SysSolver.m

clear, close all; clc

tspan = [0, 10];
initConds = [pi/2, 0];
timeStep = 0.1;
[t, y] = rk4sys(@dydtSys, tspan, initConds, timeStep);

y1 = y(:, 1);
y2 = y(:, 2);

figure
plot(t, y(:, 1), '-o', t, y(:, 2), '-s')
xlabel('\bf t')
legend('y_{1}(t)', 'y_{2}(t)')
grid on