# Need help solving a system of 2 second-order differential equations using ode45

43 views (last 30 days)
Caleb Rudick on 22 Feb 2021
Edited: Star Strider on 22 Feb 2021
I need help fixing my code to solve the two second-order equations.
y'' = gcos(theta) + (Ln+y)(theta')^2 - ky/m;
theta'' = (-gsin(theta) - 2(theta')(y'))/(Ln+y);
I can't find what's wrong with my code but I'm not very familiar with ode45 and my plot doesn't look right. What I have to do is create a plot of theta over time.
This is what I should be getting:
This is what I've been getting:
Here's my code:
syms y(t) theta(t)
% Define constants
g = 9.8; % m/s^2
m = 1; % kg
k = 100; % N/m
Ln = .5; % m
% Define initial conditions
y0 = .2;
ydot0 = 0;
theta0 = 1; % degrees
initialVals = [.2,1,0,0];
% Give equations
eqn1 = diff(y,t,2) == g*cosd(theta) + (Ln+y)*diff(theta,t)^2 - k*y/m;
eqn2 = diff(theta,t,2) == (-g*sind(theta)-2*diff(theta,t)*diff(y,t))/(Ln+y);
% Solve the system
vars = [y(t);theta(t)];
eqnSystem = odeToVectorField([eqn1,eqn2]);
eqnSystem = matlabFunction(eqnSystem,'vars', {'t','Y'});
thetaVals = ode45(eqnSystem,[0,16],initialVals);
tVec = linspace(0,16);
thetaVec = deval(thetaVals,tVec,1);
plot(tVec,thetaVec);
title('When Ln = .5m')
xlabel('Time (seconds')
ylabel('Pendulum angle (degrees)')
Star Strider on 22 Feb 2021
What about it does not work as you want it to?
When I ran it, it ran without error and produced the appropriate plot.

Star Strider on 22 Feb 2021
Edited: Star Strider on 22 Feb 2021
I am not certain what you want.
This assignment:
thetaVec = deval(thetaVals,tVec,1);
returns the value of ‘theta’ for the times defined by ‘tVec’, and then:
plot(tVec,thetaVec)
plots it appropriately.
If you want to plot all the solutions:
[t,y] = ode45(eqnSystem,[0,16],initialVals);
figure
plot(t, y);
title('When Ln = .5m')
xlabel('Time (seconds')
ylabel('Pendulum angle (degrees)')
legend(string(Subst), 'Location','SW')
will plot all of them with legend entries for all of them.
EDIT — (22 Feb 2021 at 18:30)
One problem is that ‘initVals’ is incorrect. It should be:
initialVals = [theta0,thetadot0,y0,ydot0]; % Use the ‘Subs’ Output From ‘odeToVectorField’ To Determine How It Assigned These
Unfortunately, none of the outputs of your differential equation system look like what you’re supposed to get, even with that correction. I don’t see any errors in the way you’ve coded the differential equation systems themselves. Check to be certain that the constants are correct.