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

6 次查看(过去 30 天)
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

回答(2 个)

Angelo Yeo
Angelo Yeo 2023-7-20
编辑:Angelo Yeo 2023-7-20
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.

David Goodmanson
David Goodmanson 2023-7-20
编辑:David Goodmanson 2023-7-20
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

类别

Help CenterFile Exchange 中查找有关 Argument Definitions 的更多信息

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by