MATLAB Answers

Please help i seem to be failing to understand what the following errors mean and how to correct for them.

2 views (last 30 days)
yolanda dube
yolanda dube on 2 Sep 2021
Commented: Matt J on 2 Sep 2021
%Script file for the Analysis of the double Pendulum(Chaos)
%Purpose : School Assignment
%written by : Yolanda Dube
g = 9.8; %gravitational Acceleration
l1 = 1; %length of the first rope
l2 = 1; %lenght of the second rope
m1 = 0.5; %mass of particle 1
m2 = 0.5; %mass of particle 2
t = 0:0.02:20*pi; %time period
th1 = 45; %angular displacement
th2 = 0; %angular displacement
v1 = 0; %initial velocity of the first particle
v2 = 0; %initial velocity of the second particle
yi = (pi/180).*[th1;0;th2;0];
%%using the differential equation tool
t0 = 0;
[t,y] = ode45(@yi,t0,t,f)
%%trajectory
x1 = l1*sin(y(1,:));
y1 = l1*cos(y(1,:));
x2 = l1*sin(y(1,:)) + 12*sin(y(3,:));
y2 = -l1*cos(y(1,:)) - 12*cos(y(3,:));
subplot(2,2,1)
xgrid
plot(x1,y1,'-ob')
plot(x2,y2,'-dr')
xtitle('TRAJECTORY')
subplot(2,2,2)
xgrid
plot(t, (180/pi)*y(1,:),'-ob')
plot(t,(180/pi)*y(3,:),'-dr')
xtitle('Angular displacement of m1 AND m2')
subplot(2,2,3)
xgrid
plot(t,(180/pi)*y(2,:),'-ob')
plot(t,(180/pi)*y(4,:),'-dr')
xtitle('Velocity of m1 AND m2')
subplot(2,2,4)
xgrid
plot(t,(180/pi).*y(1,:),(180/pi).*y(3,:),'k')
%%Looping our sytems of equations
function [ydash]=f(t,y)
ydash(t,1)=y(2)
ydash(t,2)=(1/((m1+m2)*l1-m2*l1*((cos(y(1)-y(3)))^2)))*(-m2*l1*(y(2)*y(2))*sin(y(1)-y(3))*cos(y(1)-y(3)) + m2*g*sin(y(3))*cos(y(1)-y(3)) - m1*l2*(y(4)*y(4))*sin(y(1)-y(3)) - (m1+m2)*g*sin(y1))
ydash(t,3) = y(4)
y(4) = (1/(l2-(m2*l2/(m1+m2))*((cos(y(1)-y(3)))^2)))*(l1*(y(2)*y(2))*sin(y(1)-y(3)) - g*sin(y(3)) + (m2*l2/(m1+m2))*(y(4)*y(4))*sin(y(1)-y(3))*cos(y(1)-y(3)) + g*sin(y(1))*cos(y(1)-y(3)))
end
  4 Comments
Stephen
Stephen on 2 Sep 2021
yi is a 4x1 numeric, you cannot define a function handle from it by simply placing '@' in front of its variable name.
How to define valid function handles is covered quite well in the MATLAB documentation:
The fourth input argument to ODE45 is defined in the documentation as a structure of options for the ODE solver. Instead you called the function f without any input arguments (thus the error message you are getting). Even if you provided f with input arguments it does not return a structure of options, so it is not clear why you expect this to be a valid fourth input to ODE45.
I suspect that you should be providing a handle to f as the first input.
It is not clear what you expect defining the value of y at the very end of f is going to achieve.

Sign in to comment.

Accepted Answer

Stephen
Stephen on 2 Sep 2021
Edited: Stephen on 2 Sep 2021
I suspect that something like this is what you are trying to do (I fixed many small bugs):
g = 9.8; %gravitational Acceleration
l1 = 1; %length of the first rope
l2 = 1; %lenght of the second rope
m1 = 0.5; %mass of particle 1
m2 = 0.5; %mass of particle 2
th1 = 45; %angular displacement
th2 = 0; %angular displacement
v1 = 0; %initial velocity of the first particle
v2 = 0; %initial velocity of the second particle
y0 = (pi/180).*[th1;0;th2;0];
ts = 0:0.02:20*pi; %time period
fun = @(t,y)f(t,y,m1,m2,l1,l2,g); % parameterize function
[t,y] = ode45(fun,ts,y0)
t = 3142×1
0 0.0200 0.0400 0.0600 0.0800 0.1000 0.1200 0.1400 0.1600 0.1800
y = 3142×4
0.7854 0 0 0 0.7835 -0.1849 0.0013 0.1310 0.7780 -0.3704 0.0052 0.2636 0.7687 -0.5572 0.0119 0.3997 0.7557 -0.7460 0.0213 0.5409 0.7389 -0.9371 0.0335 0.6884 0.7182 -1.1312 0.0488 0.8438 0.6936 -1.3295 0.0673 1.0089 0.6650 -1.5318 0.0893 1.1839 0.6323 -1.7372 0.1148 1.3673
function ydash = f(t,y,m1,m2,l1,l2,g)
ydash = nan(4,1); % preallocate!
ydash(1) = y(2);
ydash(2) = (1/((m1+m2)*l1-m2*l1*((cos(y(1)-y(3)))^2)))*(-m2*l1*(y(2)*y(2))*sin(y(1)-y(3))*cos(y(1)-y(3)) + m2*g*sin(y(3))*cos(y(1)-y(3)) - m1*l2*(y(4)*y(4))*sin(y(1)-y(3)) - (m1+m2)*g*sin(y(1)));
ydash(3) = y(4);
ydash(4) = (1/(l2-(m2*l2/(m1+m2))*((cos(y(1)-y(3)))^2)))*(l1*(y(2)*y(2))*sin(y(1)-y(3)) - g*sin(y(3)) + (m2*l2/(m1+m2))*(y(4)*y(4))*sin(y(1)-y(3))*cos(y(1)-y(3)) + g*sin(y(1))*cos(y(1)-y(3)));
end
See also:

More Answers (1)

Matt J
Matt J on 2 Sep 2021
Edited: Matt J on 2 Sep 2021
Not sure what you are trying to do with ode45. The way you are invoking it does not match any syntax in the documentation. Perhaps this is what you meant,
[t,y] = ode45(@f,t,yi)

Tags

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!