Question

1. Use Matlab to solve the differential equation (d^2φ/dt)=-(g/R)sin(φ), for the case that the board is released from φ0 = 20 degrees, using the values R = 5 m and g = 9.8 m/s^2 . Make a plot of φ against time for two or three periods. To do this, you'll need two .m files: one with your main code, which calls ode45, and one with the differential equation you're solving.

2. On the same picture, plot the approximate solution φ(t)=φ0cos(⍵t), with the same φ0 = 20 degrees. (⍵=sqrt(g/R))

3. Generate two new plots comparing the numerical and analytical solutions for φ0 = 40 degrees and φ0 = 90 degrees.

Answer 1

clc;clear all;close all;
%initial condition phi=20degree
tspan = [0:0.01:50];
y0=[20;0];
[t, y] = ode45('fun_1', tspan, y0);
omega=sqrt(9.8/5);
y1=20*cos(omega*tspan);
figure (1)
plot(t,y(:,1));
hold on
plot(t,y1);
xlabel('t');
ylabel('phi(t)in degree');
title('initial conditon phi=20degree')
%initial condition phi=40degree
y0=[40;0];
[t, y] = ode45('fun_1', tspan, y0);
omega=sqrt(9.8/5);
y2=40*cos(omega*tspan);
figure (2)
plot(t,y(:,1));
hold on
plot(t,y2);
xlabel('t');
ylabel('phi(t)in degree');
title('initial conditon phi=40degree')
%initial condition phi=90degree
y0=[90;0];
[t, y] = ode45('fun_1', tspan, y0);
omega=sqrt(9.8/5);
y3=90*cos(omega*tspan);
figure (3)
plot(t,y(:,1));
hold on
plot(t,y3);
xlabel('t');
ylabel('phi(t)in degree');
title('initial conditon phi=90degree')

%the file name should be in fun_1
function f = fun_1(t,y)
f=zeros(2,1);
%state space form of the differential equation
f(1)=y(2);
f(2)=-(9.8/5)*sind(y(1));