Trying to use modified secant method to solve for a theta over different g values.

When I run the code below I am left with NaN for all the final_theta values. Before this I was getting answers that were in the 10^-321 which is also wrong. Any explanation for this woul be helpful.

clear
clc
close all
func = @(d, theta, g, v, y) tan(theta) * d - (g / (2 * v^2 * (cos(theta))^2)) * d^2 + y; 
y = 0; %m
d = 90; % m
v = 30; % m/s 
d_theta = .01; %
theta = 1.38; %initial guess for theta in degrees between [0,90]
if theta > 1.58 || theta < 0
    error("invalid initial guess")
end
epf = .01; %desired error
celestial_bodies = ["Mercury", "Venus", "Earth", "Moon", "Mars", "Saturn", "Uranus"];
relative_g = [0.378, 0.907, 1.000, 0.166, 0.377, 0.916, 0.889];
final_theta = zeros(length(celestial_bodies));
for i = 1:length(celestial_bodies)
    g = relative_g(i) * 9.81;
    [final_theta(i), ~, ~] = modified_secant(func, theta, epf, d_theta, d, g, y, v);
    final_theta(i) = rad2deg(final_theta(i));
    disp([celestial_bodies(i) ' Angle theta :' num2str(final_theta(i)) 'degrees'])
end
function [root, ep, n] = modified_secant(func, xs, epf, dx, varargin)
% dx : fractional perturbation of the sole initial guess
n = 0;
ep = 100;
while ep > epf
    xs_new = xs - ( dx * func(xs, varargin{:})) / ( func(xs, varargin{:}) - func((xs - dx), varargin{:}));
    n = n + 1;
    if xs_new ~= 0
        ep = abs( (xs_new - xs) / xs_new) * 100;
    end
    xs = xs_new;
end
root = xs_new;
end

Answer 1

[Angelo Yeo]

You can debug the script. See that the input v for func in the first iteration is equal to 0. This makes the output from the func in the first iteration -Inf.

If you are not familiar with how to debug, check out the documentation below.

https://www.mathworks.com/help/matlab/matlab_prog/debugging-process-and-features.html

Answer 2

[David Goodmanson]

Hi Carson,

let b = g*d^2/(2*v^2)   % a distance

then your equation (multiplied by -1) is

b/cos^2 - d*tan - y = 0
plug in    1/cos^2 = (1 + tan^2)   then
b*(1 + tan^2) - d*tan -y = 0

which is a quadratic in tan. you can solve for that.

d = 90; % m
v = 30; % m/s 
g = 9.81;
b = g*d^2/(2*v^2);
tanth = roots([b,-d,b-y])
theta = atand(tanth)          % degrees
   theta =
     54.6496
     32.1706
% check original equation for both values of theta
check = tand(theta)*d - (g*d^2./(2*v^2*(cosd(theta)).^2)) + y*ones(2,1)
  check =
     1.0e-13 *
           0
      0.1421