error that keeps happening
조회 수: 1 (최근 30일)
이전 댓글 표시
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')
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');
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
댓글 수: 1
Walter Roberson
2023년 5월 3일
clear all; close all;
in the middle of your code is not doing you any favors.
답변 (1개)
Image Analyst
2023년 5월 4일
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
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!