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

4 次查看（过去 30 天）
Carson
Carson 2023-7-19
编辑： David Goodmanson 2023-7-20
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.
https://www.mathworks.com/help/matlab/matlab_prog/debugging-process-and-features.html
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 Center 和 File Exchange 中查找有关 Applications 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by