Week 3 - Solving second order ODEs

Function script

function [dtheta_dt] = ode_fun(t,theta,b,g,l,m)

theta1 = theta(1)
theta2 = theta(2)
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2 - (g/l)*sind(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

main script

clear 
close
clc

%defining the inputs to variable
b=0.05;
g=9.81;
l=0.75;
m=1;
%initial ciondition angular displacement and angular velocity
theta_0=[0;5];
%time inputs
tspan=linspace(0,10,200);

%calling the cunnction of ode
[t,result]=ode45(@(t,theta) ode_fun(t,theta,b,g,l,m),tspan,theta_0);

Theta1 = result(:,1);
ct=1;
%pendulum animation 
for i = 1:length(Theta1)
    THETA1 = Theta1(i);
   
        
x0=0;
y0=0;

x1 = l*sin(THETA1);
y1 = -l*cos(THETA1);


%ploting

plot([x0 x1],[y0 y1],'linewidth',3)

hold on
plot(x1,y1,'o','markerfacecolor',[1 0.5 0],'MarkerSize',20)
axis([-1.5 1.5 -1.5 1.5])
M(ct) = getframe(gcf);
ct=ct+1;
hold off

end

movie(M)
videofile = VideoWriter('pendulum_od.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Steps used in solving the ODE problem:

Step1. Firstly we discretise the ode function on rough paper , which makes 2 part of equation 1^st one is dtheta1_dt , and 2^nd one is dtheta2_dt. Which is look like the equation given below :

dtheta1_dt = theta2;

dtheta2_dt = -(b/m)*theta2 - (g/l)*sind(theta1);

step2. Moving to the next step we defines the parameters like B,G,L,M,

step3.then we create the theta array theta_0

step4. Next we gives the ‘ linspace ‘ command to tspan or timespan in which first is velocity and second is displacement and last is number of values.

tspan=linspace(0,10,200);

step5.

1. Now its time to create the function “ode_fun” , opening the new script we gives the function command to create our above function which we gives in first step. In function We provides the output and input first which is lokes like this [dtheta_dt] = ode_fun(t,theta,b,g,l,m)
2. after giving inputs we defines them below

3.now our discetise equation comes in play we have to give thatb equation here and later arrange that equation like wise given in question

dtheta1_dt = theta2;

dtheta2_dt = -(b/m)*theta2 - (g/l)*sind(theta1);

dtheta_dt = [dtheta1_dt;dtheta2_dt];

4. now end the function here and save it.

Step6. Now here is our most important step which is calling the function

[t,result]=ode45(@(t,theta) ode_fun(t,theta,b,g,l,m),tspan,theta_0);

1.t is the column vector If tspan contains two elements [t0 tf], then t contains the internal evaluation points used to perform the integration.

2.result is column vector

3.ode45 is a versatile ODE solver and is the first solver you should try for most problems. However, if the problem is stiff or requires high accuracy, then there are other ODE solvers that might be better suited to the problem.

4. @ is use to call back the function . @ is provided before the function.Indirectly calling a function enables you to invoke the function regardless of where you call it from.
5. ode_fun(t,theta,b,g,l,m) is our function

Step7. Now next is creating the pendulum animation using the for loop we can create it in for loop we gives the code given below , in which we use theta1 which we defines outside the for loop and as we creating it for the moving action hence we need the multiple values for which we use our results 1^st column which is our velocity column

Theta1 = result(:,1);

For i = 1:length(Theta1)

THETA1 = Theta1(i);

x0=0;

y0=0;

x1 = l*sin(THETA1);

y1 = -l*cos(THETA1);

%ploting

plot([x0 x1],[y0 y1],'linewidth',3)

hold on

plot(x1,y1,'o','markerfacecolor',[1 0.5 0],'MarkerSize',20)

axis([-1.5 1.5 -1.5 1.5])

M(ct) = getframe(gcf);

ct=ct+1;

hold off

end

in above for loop we next defines the inputs x1 y1 and x0 y0 . as our one point of pendulum is fixed . wer gives x0=o; and y0=0; . as our one point is moving we gives the

x1 = 1*sin(theta1)

y1 = 1*cos(theta1)

step8. next step we use here is “ ploting command “ 1^st command is to plot the line or penulum string and second plot command is to make the ball of pendulum by using “ marker command “ In between first plot and second plot we use the hold on command to hold both plot till end later we make it off also giving the hold off command

step9. Next command is to get the frames for required movie file for that we gives the command m(ct)=getframe(gcf);

step10. Our last and final command is to make the movie of the created pendulum for which we gives the setn of command given below . in which first we create videofile and then give a name to movie file

movie(M)

videofile = VideoWriter('pendulum_od.avi','Uncompressed AVI');

open(videofile)

writeVideo(videofile,M)

close(videofile)

ERROR FACED

"conversion to double fdrom struct is not possible" this error occur because we give same variable diffrent time for diffrent commands . first we use m as mass and second time we use it for frame creating command which is m(ct). so to enhanced tis error we make a slight chage in which we make smal m to capital M in m(ct) command which makes it M(CT). and our error solved

animation video file

https://drive.google.com/file/d/1qrricAQrWpPV_gJEkqaCJth8XfBOPmUD/view?usp=sharing

note: i am unable to upload the video file here facing some issues hence i uploaded it to google drive to check the animation video please refer above link