error that keeps happening

When I run this code everything works up until Part d where this error keeps showing up, although I did define the function earlier:

Unrecognized function or variable 'theta_explicit'.

Error in hw4_405758924_p1 (line 138)

plot(t,theta_explicit,'b',t,theta_semi_implicit,'r',t,theta_implicit,'g');

%% Part a

% Pendulum physics problem using explicit Euler method
clear all; close all;

% Define constants
g = 9.8;        % gravitational acceleration (m/s^2)
L = 2.5;        % length of string (m)
dt = 0.01;      % time step (s)
t_end = 20;     % end time (s)

% Initial conditions
theta0 = pi/2;  % initial angle (radians)
w0 = 0;         % initial angular velocity (radians/s)

% Pre-allocate arrays
N = round(t_end/dt);
theta_explicit = zeros(N+1,1);
w_explicit = zeros(N+1,1);
E_total = zeros(N+1,1);

% Set initial conditions
theta_explicit(1) = theta0;
w_explicit(1) = w0;
E_total(1) = g*L*(1-cos(theta_explicit(1)));

% Compute theta and w at each time step
for k = 1:N
    w_explicit(k+1) = w_explicit(k) - g*L*sin(theta_explicit(k))*dt;
    theta_explicit(k+1) = theta_explicit(k) + w_explicit(k)*dt;
    E_total(k+1) = g*L*(1-cos(theta_explicit(k+1))) + 0.5*L^2*w_explicit(k+1)^2;
end

% Plot results
t = linspace(0,t_end,N+1);
plot(t,theta_explicit,'b',t,w_explicit,'r',t,E_total,'g');
title('Explicit Euler Method');
xlabel('Time (s)');
ylabel('Angular Position (radians) / Angular Velocity (radians/s) / Total Energy per unit mass (J/kg)');
legend('Angular Position','Angular Velocity','Total Energy');
saveas(gcf,'explicit_euler.png');

%% Part b

% Pendulum physics problem using semi-implicit Euler method
clear all; close all;

% Define constants
g = 9.8;        % gravitational acceleration (m/s^2)
L = 2.5;        % length of string (m)
dt = 0.01;      % time step (s)
t_end = 20;     % end time (s)

% Initial conditions
theta0 = pi/2;  % initial angle (radians)
w0 = 0;         % initial angular velocity (radians/s)

% Pre-allocate arrays
N = round(t_end/dt);
theta_semi_implicit = zeros(N+1,1);
w_semi_implicit = zeros(N+1,1);
E_total = zeros(N+1,1);

% Set initial conditions
theta_semi_implicit(1) = theta0;
w_semi_implicit(1) = w0;
E_total(1) = g*L*(1-cos(theta_semi_implicit(1)));

% Compute theta and w at each time step
for k = 1:N
    w_semi_implicit(k+1) = w_semi_implicit(k) - g*L*sin(theta_semi_implicit(k+1))*dt;
    theta_semi_implicit(k+1) = theta_semi_implicit(k) + w_semi_implicit(k+1)*dt;
    E_total(k+1) = g*L*(1-cos(theta_semi_implicit(k+1))) + 0.5*L^2*w_semi_implicit(k+1)^2;
end

% Plot results
t = linspace(0,t_end,N+1);
plot(t,theta_semi_implicit,'b',t,w_semi_implicit,'r',t,E_total,'g');
title('Semi-Implicit Euler Method');
xlabel('Time (s)');
ylabel('Angular Position (radians) / Angular Velocity')

%% Part c
% Constants
g = 9.8; % m/s^2
L = 2.5; % m
dt = 0.01; % s
t_start = 0; % s
t_end = 20; % s
theta0 = pi/2; % rad
n_steps = ceil((t_end - t_start)/dt);
% Initialize arrays
theta_implicit = zeros(n_steps, 1);
w_implicit = zeros(n_steps, 1);
energy = zeros(n_steps, 1);
% Initial conditions
theta_implicit(1) = 5; %<<< made up value for theta0;
w_implicit(1) = 0; % starting from rest
g = 9.80665; %<<< added missing value for g (m/sec^2)
L = 1; %<<< made up value for L
% Solve using implicit Euler method
for i = 2:n_steps
    w_implicit(i) = w_implicit(i-1) - g*L*sin(theta_implicit(i-1))*dt;
    theta_implicit(i) = newton_raphson(theta_implicit(i-1), w_implicit(i), dt);
    energy(i) = g*L*(1 - cos(theta_implicit(i))) + 0.5*(L*w_implicit(i))^2;
end
% Plotting
t = linspace(t_start, t_end, n_steps);
hold on
plot(t, theta_implicit, 'r', 'LineWidth', 2)
title('Angular position of pendulum')
xlabel('Time (s)')
ylabel('Angular position (rad)')
legend('Explicit', 'Semi-implicit', 'Implicit')
saveas(gcf, 'pendulum_position.png')

figure()
hold on
plot(t, energy, 'r', 'LineWidth', 2)
title('Total energy of pendulum')
xlabel('Time (s)')
ylabel('Energy per unit mass (J/kg)')
legend('Explicit', 'Semi-implicit', 'Implicit')

Warning: Ignoring extra legend entries.

saveas(gcf, 'pendulum_energy.png')

%% Part d 

% Plot angular position
figure;
plot(t,theta_explicit,'b',t,theta_semi_implicit,'r',t,theta_implicit,'g');

Unrecognized function or variable 'theta_explicit'.

title('Angular Position of Mass vs Time');
xlabel('Time (s)');
ylabel('Angular Position (radians)');
legend('Explicit Euler','Semi-Implicit Euler','Implicit Euler');
saveas(gcf,'angular_position.png');




%% Part e
% Calculating total energy per unit mass for each method
g = 9.8;        % gravitational acceleration (m/s^2)
E_explicit = g*L*(1 - cos(theta_explicit)) + 0.5*L^2*w_explicit.^2;
E_semi_implicit = g*L*(1 - cos(theta_semi_implicit)) + 0.5*L^2*w_semi_implicit.^2;
E_implicit = g*L*(1 - cos(theta_implicit)) + 0.5*L^2*w_implicit.^2;

% Plotting total energy per unit mass for each method
figure;
plot(t, E_explicit, 'r-', 'LineWidth', 1.5);
hold on;
plot(t, E_semi_implicit, 'b-', 'LineWidth', 1.5);
plot(t, E_implicit, 'g-', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Total Energy per unit mass (J/kg)');
title('Total Energy per unit mass of Simple Pendulum');
legend('Explicit Euler', 'Semi-Implicit Euler', 'Implicit Euler');
grid on;

% Saving the plot
saveas(gcf, 'pendulum_energy.png');

% Newton-Raphson method to solve implicit Euler equation
function [theta_next] = newton_raphson(theta_curr, omega_curr, dt)
g = 9.80665; %<<< added missing value for g (m/sec^2)
L = 1; %<<< made up value for L
theta_next = theta_curr; % initial guess
f = @(theta_next) theta_next - theta_curr - dt*omega_curr - (dt^2/2)*(-g*L*sin(theta_next));
df = @(theta_next) 1 - (dt^2/2)*(-g*L*cos(theta_next));
for i = 1:5 % 5 iterations for convergence
    theta_next = theta_next - f(theta_next)/df(theta_next);
end
end

[Walter Roberson]

clear all; close all;

in the middle of your code is not doing you any favors.

Answer 1

[Image Analyst]

In line 46 you have a clear that erases the variable:

%% Part b
% Pendulum physics problem using semi-implicit Euler method
clear all; close all;

Get rid of that clear all close all line. There were quite a few more errors in there. I believe I've fixed them for you.

%% Part a
% Pendulum physics problem using explicit Euler method
clear all; 
close all;
clc;
% Define constants
g = 9.8;        % gravitational acceleration (m/s^2)
L = 2.5;        % length of string (m)
dt = 0.01;      % time step (s)
t_end = 20;     % end time (s)
% Initial conditions
theta0 = pi/2;  % initial angle (radians)
w0 = 0;         % initial angular velocity (radians/s)
% Pre-allocate arrays
N = round(t_end/dt);
theta_explicit = zeros(N+1,1);
w_explicit = zeros(N+1,1);
E_total = zeros(N+1,1);
% Set initial conditions
theta_explicit(1) = theta0;
w_explicit(1) = w0;
E_total(1) = g*L*(1-cos(theta_explicit(1)));
% Compute theta and w at each time step
for k = 1:N-1
    w_explicit(k+1) = w_explicit(k) - g*L*sin(theta_explicit(k))*dt;
    theta_explicit(k+1) = theta_explicit(k) + w_explicit(k)*dt;
    E_total(k+1) = g*L*(1-cos(theta_explicit(k+1))) + 0.5*L^2*w_explicit(k+1)^2;
end
% Plot results
t = linspace(0,t_end,N+1)';
figure('Name', 'Part a', 'NumberTitle', 'off')
plot(t,theta_explicit,'b',t,w_explicit,'r',t,E_total,'g');
grid on;
title('Explicit Euler Method');
xlabel('Time (s)');
ylabel('Angular Position (radians) / Angular Velocity (radians/s) / Total Energy per unit mass (J/kg)');
legend('Angular Position','Angular Velocity','Total Energy');
drawnow; % Force immediate screen update.
saveas(gcf,'explicit_euler.png');
%% Part b
% Pendulum physics problem using semi-implicit Euler method
% Define constants
g = 9.8;        % gravitational acceleration (m/s^2)
L = 2.5;        % length of string (m)
dt = 0.01;      % time step (s)
t_end = 20;     % end time (s)
% Initial conditions
theta0 = pi/2;  % initial angle (radians)
w0 = 0;         % initial angular velocity (radians/s)
% Pre-allocate arrays
N = round(t_end/dt);
theta_semi_implicit = zeros(N+1,1);
w_semi_implicit = zeros(N+1,1);
E_total = zeros(N+1,1);
% Set initial conditions
theta_semi_implicit(1) = theta0;
w_semi_implicit(1) = w0;
E_total(1) = g*L*(1-cos(theta_semi_implicit(1)));
% Compute theta and w at each time step
for k = 1:N-1
    w_semi_implicit(k+1) = w_semi_implicit(k) - g*L*sin(theta_semi_implicit(k+1))*dt;
    theta_semi_implicit(k+1) = theta_semi_implicit(k) + w_semi_implicit(k+1)*dt;
    E_total(k+1) = g*L*(1-cos(theta_semi_implicit(k+1))) + 0.5*L^2*w_semi_implicit(k+1)^2;
end
% Plot results
t = linspace(0,t_end,N+1)';
figure('Name', 'Part b', 'NumberTitle', 'off')
plot(t,theta_semi_implicit,'b',t,w_semi_implicit,'r',t,E_total,'g');
grid on;
title('Semi-Implicit Euler Method');
xlabel('Time (s)');
ylabel('Angular Position (radians) / Angular Velocity')
drawnow; % Force immediate screen update.
%% Part c
% Constants
g = 9.8; % m/s^2
L = 2.5; % m
dt = 0.01; % s
t_start = 0; % s
t_end = 20; % s
theta0 = pi/2; % rad
n_steps = N+1; %ceil((t_end - t_start)/dt);
% Initialize arrays
theta_implicit = zeros(n_steps, 1);
w_implicit = zeros(n_steps, 1);
energy = zeros(n_steps, 1);
% Initial conditions
theta_implicit(1) = 5; %<<< made up value for theta0;
w_implicit(1) = 0; % starting from rest
g = 9.80665; %<<< added missing value for g (m/sec^2)
L = 1; %<<< made up value for L
% Solve using implicit Euler method
for k = 2:n_steps
    w_implicit(k) = w_implicit(k-1) - g*L*sin(theta_implicit(k-1))*dt;
    theta_implicit(k) = newton_raphson(theta_implicit(k-1), w_implicit(k), dt);
    energy(k) = g*L*(1 - cos(theta_implicit(k))) + 0.5*(L*w_implicit(k))^2;
end
% Plotting
figure('Name', 'Part c', 'NumberTitle', 'off')
t = linspace(t_start, t_end, n_steps);
hold on
plot(t, theta_implicit, 'r', 'LineWidth', 2)
grid on;
title('Angular position of pendulum')
xlabel('Time (s)')
ylabel('Angular position (rad)')
saveas(gcf, 'pendulum_position.png')
drawnow; % Force immediate screen update.
hold on
plot(t, energy, 'r', 'LineWidth', 2)
grid on;
title('Total energy of pendulum')
xlabel('Time (s)')
ylabel('Energy per unit mass (J/kg)')
legend('theta_implicit', 'Energy')
saveas(gcf, 'pendulum_energy.png')
drawnow; % Force immediate screen update.
%% Part d
% Plot angular position
figure('Name', 'Part d', 'NumberTitle', 'off')
plot(t,theta_explicit,'b',t,theta_semi_implicit,'r',t,theta_implicit,'g');
grid on;
title('Angular Position of Mass vs Time');
xlabel('Time (s)');
ylabel('Angular Position (radians)');
legend('Explicit Euler','Semi-Implicit Euler','Implicit Euler');
saveas(gcf,'angular_position.png');
drawnow; % Force immediate screen update.
%% Part e
% Calculating total energy per unit mass for each method
g = 9.8;        % gravitational acceleration (m/s^2)
E_explicit = g*L*(1 - cos(theta_explicit)) + 0.5*L^2*w_explicit.^2;
E_semi_implicit = g*L*(1 - cos(theta_semi_implicit)) + 0.5*L^2*w_semi_implicit.^2;
E_implicit = g*L*(1 - cos(theta_implicit)) + 0.5*L^2*w_implicit.^2;
% Plotting total energy per unit mass for each method
figure('Name', 'Part e', 'NumberTitle', 'off')
plot(t, E_explicit, 'r-', 'LineWidth', 1.5);
grid on;
hold on;
plot(t, E_semi_implicit, 'b-', 'LineWidth', 1.5);
plot(t, E_implicit, 'g-', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Total Energy per unit mass (J/kg)');
title('Total Energy per unit mass of Simple Pendulum');
legend('Explicit Euler', 'Semi-Implicit Euler', 'Implicit Euler');
grid on;
drawnow; % Force immediate screen update.
% Saving the plot
saveas(gcf, 'pendulum_energy.png');
uiwait(helpdlg('Done!'))

%================================================================================
% Newton-Raphson method to solve implicit Euler equation
function [theta_next] = newton_raphson(theta_curr, omega_curr, dt)
g = 9.80665; %<<< added missing value for g (m/sec^2)
L = 1; %<<< made up value for L
theta_next = theta_curr; % initial guess
f = @(theta_next) theta_next - theta_curr - dt*omega_curr - (dt^2/2)*(-g*L*sin(theta_next));
df = @(theta_next) 1 - (dt^2/2)*(-g*L*cos(theta_next));
for k = 1:5 % 5 iterations for convergence
    theta_next = theta_next - f(theta_next)/df(theta_next);
end
end