Using Euler's Method in Matlab

First time post here. Pretty frustrated right now working on this assignment for class.

Basically, the idea is to use Euler's method to simulate and graph an equation of motion. The equation of motion is in the form of an ODE.

My professor has already put down some code for slightly similar system and would like us to derive the equation of motion using Lagrange. I believe that I have derived the EOM correctly, however I am running into problems on the Matlab side of things.

What's weird is that using a similar technique on another, seperate EOM, I have no issues. So I am unsure what I am doing wrong.

Here's the code for the part that is working correctly:

close all; clear all; clc; 
% System parameters
w = 2*pi;
c = 0.02;
% Time vectors
dt = 1e-5;
t = 0:dt:4;
theta = zeros(size(t));
thetadot = zeros(size(t));

% Initial conditions
theta(1)=pi/2; %theta(0)
thetadot(1)=0; %thetadot(0)

for I = 1 : length(t)-1;
    thetaddot = -c*thetadot(I)-w^2*sin(theta(I));
    thetadot(I+1)=thetadot(I)+thetaddot*dt;
    theta(I+1)=theta(I)+thetadot(I)*dt ;
end
figure(1);
plot(t,theta,'b');
xlabel('time(s)');
ylabel('theta');
title('Figure 1');
zoom on;
% Output the plot to a pdf file, and make it 6 inches by 4 inches
printFigureToPdf('fig1.pdf', [6,4],'in');
% Open the pdf for viewing
open fig1.pdf

Everything runs fine, except Matlab complains about the printFigureToPdf command.

Now, here is the code for the problem that I am having issues with.

close all; clear all; clc; clf
% System parameters
m=0.2;
g=9.81;
c=.2;
d=0.075;
L=0.001; %L is used for Gamma
B=0.001; %B is used for Beta
W=210*pi; %W is used for Omega

%Time vectors
dt = 1e-6; %Time Step
t=0:dt:10; %Range of times that simulation goes through
x=zeros(size(t));
xdot=zeros(size(t));

%Initialconditions
x(1)=0;%x(0)
xdot(1)=0; %xdot(0)

for I = 1 : length(t)-1;
    xddot =-1/m*(c*xdot(I)-c*L*W*cos(W)+m*g-3*B*((d+x-L*W*sin(W*t)).^(-4)-(d-x-L*W*sin(W*t)).^(-4)));
    xdot(I+1)=xdot(I)+xddot*dt;
    x(I+1)=x(I)+xdot(I+1)*dt ;
end
figure(1);
plot(t,x,'b');
xlabel('time(s)');
ylabel('distance(m)');
title('Figure 2');
zoom on;
% Output the plot to a pdf file, and make it 6 inches by 4 inches
printFigureToPdf('fig1.pdf', [6,4],'in');
% Open the pdf for viewing
open fig1.pdf

With this code, I followed the same procedure and is giving an error on line 23: "In an assignment A(I) = B, the number of elements in B and I must be the same."

Like I said, I am confused because the other code worked okay, and this second set of code gives an error.

If anyone could give me a hand with this, I would greatly appreciate it.

Thanks in advance, Dave

Edit: As suggested by a friend, I changed

x(I+1)=x(I)+xdot(I+1)*dt to 

x(I+1)=x(I)+xdot(I)*dt

However, I am still getting an error for line 23: "In an assignment A(I) = B, the number of elements in B and I must be the same."

Line 23 is:

xdot(I+1)=xdot(I)+xddot*dt;

So, I tried adjusting the code as suggested for the other line to

xdot(I+1)=xdot(I)+xddot(I)*dt;

After making this change, Matlab gets stuck, I tried letting it run for a few minutes but won't execute. I ended up having to close and reopen the application.

Answer 1

[Geoff Hayes]

David - the problem with line

 xdot(I+1)=xdot(I)+xddot*dt;

is that the xddot is 1x10000001 vector rather than a scalar. So the code is trying to assign a 1x10000001 array to a single element in your xdot. Unlike your previous example, the xddot depends on the vector t (and x) in a couple of places

 xddot =-1/m*(c*xdot(I)-c*L*W*cos(W)+m*g-3*B*((d+x-L*W*sin(W*t)).^(-4)-(d-x-L*W*sin(W*t)).^(-4)));

When the above line is evaluated, the entire t and x vectors are used and so the xddot becomes a vector too.

I suspect that you may want to use the particular time value, t(I) and x(I) instead. So try replacing the above equation with

 xddot =-1/m*(c*xdot(I)-c*L*W*cos(W)+m*g-3*B*((d+x(I)-L*W*sin(W*t(I))).^(-4)-(d-x(I)-L*W*sin(W*t(I))).^(-4)));

and see what happens. When I ran your code with the above change, the following was produced

[David]

Thank you very much, Geoff. I will give that change a try. The equation of motion is supposed to be an oscillation of sorts, so I believe part of the problem lies with my equation that I derived with Lagrange.