Hello, I'm trying to create a plot for an equation of motion defined as a function.

I am trying to create a plot for this function but when I put 'end' before the initial conditions I get an error that reads "Function definitions in a script must appear at the end of the file. Move statements to before the function definitions." I've tried to move the initial conditions to be before the function definition and that didn't run properly. I've also tried to move 'end' to the end of the code and the plots didn't show up at all. I would appreciate any help for this. Thanks!

function dXdt=eom(~,X)

%%First Conditions

theta=X(1); %Position HC

thetadot=X(2); %Velocity HC

phi=X(3); %Position Disk

phidot=X(4); %Velocity Disk

%EOMS

thetaddot=-(g/R)*sin(theta); %Half cylinder

phiddot=-(g/r)*sin(phi); % Disk

dXdt=[thetadot;thetaddot;phidot;phiddot];

end

%Initial conditions

thetazero=0;

phizero=pi/18;

Xzero=[thetazero;0;phizero;0];

%timeframe

tspan=[0,6];

%solution

[t,X]=ode45(@eom,tspan,Xzero);

%plot

figure

subplot(2,1,1);

plot(t,X(:,1));

xlabel('Time')

ylabel('\theta')

title('Angular Position of Half Cylinder')

subplot(2,1,2);

plot(t,X(:,3));

xlabel('Time')

ylabel('\phi')

title('Angular Position of Disk')

hold on

[Les Beckham]

You haven't defined g, R, and r that are used in your eom function.

[Amanda]

I'm sorry I forgot to include that portion of code. I'm solving multiple problems and the code is long so I only posted the snippet that wasn't working. The following is the portion where I defined the variables:

% Givens

M=578; % in grams

m=162; % in grams

R=15.625; % in cm

r=3.8; % in cm

g=9.81; % m/s^2

Answer 1

[Cris LaPierre]

It sounds like you have your function inside a larger script. In your version of MATLAB, the function must be moved to the bottom of the script. You can still call the function in your script. The definition just needs to be at the bottom. Something like this.

%Initial conditions
thetazero=0;
phizero=pi/18;
Xzero=[thetazero;0;phizero;0];
%timeframe
tspan=[0,6];
%solution
[t,X]=ode45(@eom,tspan,Xzero);
%plot
figure
subplot(2,1,1);
plot(t,X(:,1));
xlabel('Time')
ylabel('\theta')
title('Angular Position of Half Cylinder')
subplot(2,1,2);
plot(t,X(:,3));
xlabel('Time')
ylabel('\phi')
title('Angular Position of Disk')
hold on

function dXdt=eom(~,X)
% Givens
M=578; % in grams
m=162; % in grams
R=15.625; % in cm
r=3.8; % in cm
g=9.81; % m/s^2

%%First Conditions
    theta=X(1); %Position HC
    thetadot=X(2); %Velocity HC
    phi=X(3); %Position Disk
    phidot=X(4); %Velocity Disk
    %EOMS
    thetaddot=-(g/R)*sin(theta); %Half cylinder
    phiddot=-(g/r)*sin(phi); % Disk
    dXdt=[thetadot;thetaddot;phidot;phiddot];
end

[Amanda]

Thank you, that removed all the errors!

[Cris LaPierre]

Here's your full code. It just needed a little clean up. Hopefully, by comparing this to the code you have, you can identify the changes I made. Note that your equations of motion do not change, just the values they are using. It is therefore not necessary to define a new function each time.

% Givens
M=578; % in grams
m=162; % in grams
R=15.625; % in cm
r=3.8; % in cm
g=9.81; % m/s^2

% HALF-CYLINDER NO DISK
%thetazero=pi/9, thetadot=0, ten seconds
%initial conditions
thetazero=pi/9; %radians
thetadotzero=0; %rad/s

%ADD COM ANALYTICAL SOLUTION
%moment of inertia
ihc=(1/4)*M*R^2;
%natural frequency
nfhc=sqrt(g/R);
%time
t=linspace(0,10,1000);
%angular position analytical
thetat=thetazero*cos(nfhc*t);

%plot
figure;
plot(t,thetat);
xlabel('Time')
ylabel('Angular Position')
title('Half-Cylinder, No Disk')
grid on

%FULL CYLINDER WITH DISK
%plot when phizero=pi/9 and phidot=0 for 10 seconds
%initial conditions
phizero=pi/9;
phidotzero=0;

%ADD COM ANALYTICAL SOLUTION
%inertia
id=(1/2)*m*r^2;
ic=(1/2)*M*R^2; 
%natural frequency
nfd=sqrt(g/r);
%Time
t=linspace(0,10,1000);
%angular displacement
phit=phizero*cos(nfd*t);

%plot
figure
plot(t,phit)
xlabel('Time')
ylabel('Angular Position')
title('Full Cylinder with Disk')
grid on

%Initial conditions
thetazero=0;
phizero=pi/18;
Xzero=[thetazero;0;phizero;0];
%timeframe
tspan=[0,6];
%solution
[t,X]=ode45(@eom,tspan,Xzero,[],M,m,R,r,g);

%plot
figure
subplot(2,1,1);
plot(t,X(:,1));
xlabel('Time')
ylabel('\theta')
title('Angular Position of Half Cylinder')
subplot(2,1,2);
plot(t,X(:,3));
xlabel('Time')
ylabel('\phi')
title('Angular Position of Disk')

%Initial conditions
thetazero=pi/18;
phizero=pi/9;
Xzero=[thetazero;0;phizero;0];
%timeframe
tspan=[0,6];
%solution
[t,X]=ode45(@eom,tspan,Xzero,[],M,m,R,r,g);

%plot
figure
subplot(2,1,1);
plot(t,X(:,1));
xlabel('Time')
ylabel('\theta')
title('Angular Position of Half Cylinder')
subplot(2,1,2);
plot(t,X(:,3));
xlabel('Time')
ylabel('\phi')
title('Angular Position of Disk')

%Initial conditions
thetazero=pi/9;
phizero=-pi/9;
Xzero=[thetazero;0;phizero;0];
%timeframe
tspan=[0,6];
%solution
[t,X]=ode45(@eom,tspan,Xzero,[],M,m,R,r,g);

%plot
figure
subplot(2,1,1);
plot(t,X(:,1));
xlabel('Time')
ylabel('\theta')
title('Angular Position of Half Cylinder')
subplot(2,1,2);
plot(t,X(:,3));
xlabel('Time')
ylabel('\phi')
title('Angular Position of Disk')

function dXdt=eom(~,X,M,m,R,r,g)
%%First Conditions
    theta=X(1); %Position HC
    thetadot=X(2); %Velocity HC
    phi=X(3); %Position Disk
    phidot=X(4); %Velocity Disk
    %EOMS
    thetaddot=-(g/R)*sin(theta); %Half cylinder
    phiddot=-(g/r)*sin(phi); % Disk
    dXdt=[thetadot;thetaddot;phidot;phiddot];
end

[Amanda]

Thank you, I was able to see the changes you made and fix it!

Answer 2

[Steven Lord]

Please show (or even attach) the full and exact file. My suspicion is that the lines where you define those constants are ahead of the function declaration, but that you forgot to actually call the function you defined, something like:

x = 1

x = 1

function y = timestwo(x)
y = 2*x;
surf(peaks)
end

While I've defined the timestwo function, I have not called that function so all this code does is define x to be 1. So there's no graphics plotted when I run this code.

If I'd put y = timestwo(x) immediately before or after the line x = 1 in the code, then it would perform the multiplication and create the peaks plot.

[Amanda]

This is the full code. Right now I'm solving five different problems and the first two portions ("Half-Cylinder No Disk" and "Full Cylinder with Disk") work fine but when I get into the functions the plots don't show up. Moving the function definition to the bottom has removed the error I was seeing, but the plots still don't show.

% Givens

M=578; % in grams

m=162; % in grams

R=15.625; % in cm

r=3.8; % in cm

g=9.81; % m/s^2

% HALF-CYLINDER NO DISK

%thetazero=pi/9, thetadot=0, ten seconds

%initial conditions

thetazero=pi/9; %radians

thetadotzero=0; %rad/s

%ADD COM ANALYTICAL SOLUTION

%moment of inertia

ihc=(1/4)*M*R^2;

%natural frequency

nfhc=sqrt(g/R);

%time

t=linspace(0,10,1000);

%angular position analytical

thetat=thetazero*cos(nfhc*t);

%plot

figure;

plot(t,thetat);

xlabel('Time')

ylabel('Angular Position')

title('Half-Cylinder, No Disk')

grid on

hold on

%FULL CYLINDER WITH DISK

%plot when phizero=pi/9 and phidot=0 for 10 seconds

%initial conditions

phizero=pi/9;

phidotzero=0;

%ADD COM ANALYTICAL SOLUTION

%inertia

id=(1/2)*m*r^2;

ic=(1/2)*M*R^2;

%natural frequency

nfd=sqrt(g/r);

%Time

t=linspace(0,10,1000);

%angular displacement

phit=phizero*cos(nfd*t);

%plot

figure

plot(t,phit)

xlabel('Time')

ylabel('Angular Position')

title('Full Cylinder with Disk')

grid on

function dXdt=eom(~,X)

%%First Conditions

theta=X(1); %Position HC

thetadot=X(2); %Velocity HC

phi=X(3); %Position Disk

phidot=X(4); %Velocity Disk

%EOMS

thetaddot=-(g/R)*sin(theta); %Half cylinder

phiddot=-(g/r)*sin(phi); % Disk

dXdt=[thetadot;thetaddot;phidot;phiddot];

end

%Initial conditions

thetazero=0;

phizero=pi/18;

Xzero=[thetazero;0;phizero;0];

%timeframe

tspan=[0,6];

%solution

[t,X]=ode45(@eom,tspan,Xzero);

%plot

figure

subplot(2,1,1);

plot(t,X(:,1));

xlabel('Time')

ylabel('\theta')

title('Angular Position of Half Cylinder')

subplot(2,1,2);

plot(t,X(:,3));

xlabel('Time')

ylabel('\phi')

title('Angular Position of Disk')

hold on

function dXdT=eom2(~,X)

%%First Conditions

theta=X(1); %Position HC

thetadot=X(2); %Velocity HC

phi=X(3); %Position Disk

phidot=X(4); %Velocity Disk

%EOMS

thetaddot=-(g/R)*sin(theta); %Half cylinder

phiddot=-(g/r)*sin(phi); % Disk

dXdT=[thetadot;thetaddot;phidot;phiddot];

%Initial conditions

thetazero=pi/18;

phizero=pi/9;

Xzero=[thetazero;0;phizero;0];

%timeframe

tspan=[0,6];

%solution

[t,X]=ode45(@eom2,tspan,Xzero);

%plot

figure

subplot(2,1,1);

plot(t,X(:,1));

xlabel('Time')

ylabel('\theta')

title('Angular Position of Half Cylinder')

subplot(2,1,2);

plot(t,X(:,3));

xlabel('Time')

ylabel('\phi')

title('Angular Position of Disk')

hold on

end

function dXdT=eom3(~,X)

%%First Conditions

theta=X(1); %Position HC

thetadot=X(2); %Velocity HC

phi=X(3); %Position Disk

phidot=X(4); %Velocity Disk

%EOMS

thetaddot=-(g/R)*sin(theta); %Half cylinder

phiddot=-(g/r)*sin(phi); % Disk

dXdT=[thetadot;thetaddot;phidot;phiddot];

%Initial conditions

thetazero=pi/9;

phizero=-pi/9;

Xzero=[thetazero;0;phizero;0];

%timeframe

tspan=[0,6];

%solution

[t,X]=ode45(@eom3,tspan,Xzero);

%plot

figure

subplot(2,1,1);

plot(t,X(:,1));

xlabel('Time')

ylabel('\theta')

title('Angular Position of Half Cylinder')

subplot(2,1,2);

plot(t,X(:,3));

xlabel('Time')

ylabel('\phi')

title('Angular Position of Disk')

hold on

end