i keep getting this error: Error: Function definitions in a script must appear at the end of the file. Move all statements after the function definition to before the ...

i keep getting this erorr: Error: File: hw4_405758924_p1.m Line: 104 Column: 1 Function definitions in a script must appear at the end of the file. Move all statements after the "newton_raphson" function definition to before the first local function definition.

main()

Warning: Ignoring extra legend entries.

Warning: Ignoring extra legend entries.

function main
% 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 = zeros(n_steps, 1);
omega = zeros(n_steps, 1);
energy = zeros(n_steps, 1);

% Initial conditions
theta(1) = theta0;
omega(1) = 0; % starting from rest

% Solve using implicit Euler method
for i = 2:n_steps
    omega(i) = omega(i-1) - g*L*sin(theta(i-1))*dt;
    theta(i) = newton_raphson(theta(i-1), omega(i), dt);
    energy(i) = g*L*(1 - cos(theta(i))) + 0.5*(L*omega(i))^2;
end

% Plotting
t = linspace(t_start, t_end, n_steps);
hold on
plot(t, theta, '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')

% Newton-Raphson method to solve implicit Euler equation
function [theta_next] = newton_raphson(theta_curr, omega_curr, dt)
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
end

[Torsten]

Use newton_raphson as a nested function. Then your code should work (see above).

Answer 1

[Les Beckham]

You must have some code in your script after the last end that you didn't post. Make sure that the function definition is the very last thing in your script file.

Your code works after supplying values for some undefined variables.

I'm not sure why you have 3 legend entries since you are only plotting one curve in your plots. Maybe you intend to add more to the plots.

% 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 = zeros(n_steps, 1);
omega = zeros(n_steps, 1);
energy = zeros(n_steps, 1);
% Initial conditions
theta(1) = 5; %<<< made up value for theta0;
omega(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
    omega(i) = omega(i-1) - g*L*sin(theta(i-1))*dt;
    theta(i) = newton_raphson(theta(i-1), omega(i), dt);
    energy(i) = g*L*(1 - cos(theta(i))) + 0.5*(L*omega(i))^2;
end
% Plotting
t = linspace(t_start, t_end, n_steps);
hold on
plot(t, theta, 'r', 'LineWidth', 2)
title('Angular position of pendulum')
xlabel('Time (s)')
ylabel('Angular position (rad)')
legend('Explicit', 'Semi-implicit', 'Implicit')

Warning: Ignoring extra legend entries.

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')

% 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