Error plot is starting at the wrong value in my "while loop"

clear all
w=(12.35)^0.5;
%theta_exact=10*exp(-2*t)*cos(t*w)+20*exp(-2*t)*sin(w*t);
theta1(1)=10; %deg
psy1(1)=0; %deg/s
g=9.81; %m/s^2
c=4; %1/s
l=0.6; %m
dt1=0.1;
t1(1)=0;
i=1;
while t1<6  
    theta1(i+1)=theta1(i)+psy1(i)*dt1
    psy1(i+1)=(-g/l)*dt1*theta1(i)+(1-c*dt1)*psy1(i);
    theta_exact=10.*exp(-2.*t1).*cos(w.*t1)+(20/w).*exp(-2.*t1).*sin(w.*t1)
    epsilon1=(theta1(i)-theta_exact)
    t1(i+1)=t1(i)+dt1;
    i=i+1;
    
end
figure(1); clf;
hold on
plot(t1,theta1,'red')
plot(t1(1:end-1),theta_exact,'blue')
figure(2); clf;
plot(t1(1:end-1),(epsilon1))

If you read the code you notice that the second plot is a time vs "epsilon1". My issue is that the first value for "epsilon1" is -10 but epsilon1=(theta1(i)-theta_exact), where theta1(i) and theta_exact are the same value for their first values. The first value should read 0 for epsilon1, can anyone help me figure out why this is not the case?

Answer 1

[ag]

Hi @Caleb Hedin,

The issue you're encountering with the calculation of "epsilon1" arises from the order of operations within your while loop. Specifically, "epsilon1" is being calculated before "theta_exact" is updated for the next time step, which causes a mismatch in the indices used for comparison.

To resolve this, ensure that "epsilon1" is calculated after both "theta1(i)" and "theta_exact" are updated for the current time step. The below code demonstrates the required modification:

w = sqrt(12.35);
theta1(1) = 10; % deg
psy1(1) = 0; % deg/s
g = 9.81; % m/s^2
c = 4; % 1/s
l = 0.6; % m
dt1 = 0.1;
t1(1) = 0;
i = 1;

while t1(i) < 6  
    theta1(i+1) = theta1(i) + psy1(i) * dt1;
    psy1(i+1) = (-g/l) * dt1 * theta1(i) + (1 - c * dt1) * psy1(i);
    t1(i+1) = t1(i) + dt1;
    
    % Calculate theta_exact for the current time step
    theta_exact(i) = 10 * exp(-2 * t1(i)) * cos(w * t1(i)) + (20/w) * exp(-2 * t1(i)) * sin(w * t1(i));
    
    % Calculate epsilon1 for the current time step
    epsilon1(i) = theta1(i) - theta_exact(i);
    
    i = i + 1;
end

% Calculate theta_exact for the last time step
theta_exact(i) = 10 * exp(-2 * t1(i)) * cos(w * t1(i)) + (20/w) * exp(-2 * t1(i)) * sin(w * t1(i));
epsilon1(i) = theta1(i) - theta_exact(i);

% Plot results
figure(1); clf;
hold on;
plot(t1, theta1, 'red');
plot(t1, theta_exact, 'blue');
legend('Theta1', 'Theta Exact');
xlabel('Time (s)');
ylabel('Theta (degrees)');

●

figure(2); clf;
plot(t1, epsilon1);
xlabel('Time (s)');
ylabel('Epsilon (degrees)');
title('Error (Epsilon) over Time');

Hope this helps!