MATLAB Answers

Help with BVP4C solving for launch velocity with designated start and end points and variable launch angle.

3 views (last 30 days)
Hi all, I'm working on a balistic computer of sorts that will take in the initial point, the impact point, and the launch angle and will solve for the required launch velocity. Right now I'm trying to get the trejectory to match expectations without considering drag, but next steps include implementing drag and recursively finding the angle with the minimum launch velocity (which is why theta is listed as an additional parameter.
Lastly, before I get to the code, I'm fairly convinced that the problem is coming from the boundary conditions because if you look at the displayed start and end points, the start point is different than the intended start point and the angle is very different than the one that should be enforced by the fifth boundary condition. The solinit "guess" seems to have very little impact on the solution, and I'm confident with my odesys equations.
Thank you in advance for your help!
d = 100;
theta = pi/4;
solinit = bvpinit(linspace(0,d,20),@guess, theta);
sol = bvp4c(@odesys, @resfun, solinit);
fprintf('v_0 (from height) = %f \ntheta = %f \n(%f,%f) -> (%f,%f)\n', ...
sqrt(max(y(1,:))* 32.2 * 2), atand(y(2,1)/y(4,1)), y(3,1), y(1,1),y(3,end),y(1,end));
function dy = odesys(x, y, theta)
% launch parameters
g = 32.2; % gravity, ft/s^2
Cd = 0; % drag coefficient
rho = 0.00237; % density of air, slug/ft^3
c = 10; % pumpkin circumfrence (normal to motion), in
c = c / 12; % convert in to ft
r = c / 2 / pi; % radius, ft
A = pi * r ^ 2; % cross sectional area, ft^2
W = 4; % weight, lbf
m = W / g; % convert lbf to slug
% constants for use in defining dy
k = rho * Cd * A / 2 / m; % constant for drag, 1/sqrt(s*ft)
v = sqrt(y(2).^2+y(4).^2);
%defining dy
dy = y;
dy(1) = y(2);
dy(2) = -g + k * v * y(2);
dy(3) = y(4);
dy(4) = k * v * y(4);
function res = resfun(Ya, Yb, theta)
d = 100; % target distance, ft
% 50, 100, 200, and 300
h = 4; % launch height, ft
res = zeros(5,1);
% point a: (0,h)
res(1) = Ya(1)-h;
res(2) = Ya(3);
% point b: (d,0)
res(3) = Yb(1);
res(4) = Yb(3)-d;
%initial launch angle
res(5) = theta - asin(Ya(2)/sqrt(Ya(2).^2+Ya(4).^2));
function y = guess(x)
% calculate minimum launch velocity with 0 drag
g = 32.2;
sine = sind(45);
d = 100;
h = 4;
v = d / sine / sqrt( 2 * ( d + h) / g);
% guess follows 0 drag 45º launch
y(1) = -g*x^2/v^2/sine^2+x+h;
y(2) = -g*x/v/sine+v*sine;
y(3) = x;
y(4) = v*sine;
Max Friedman
Max Friedman on 20 Feb 2021
The fifth boundary condition is required to control the variable parameter. So while there are 4 DiffEqs, without the fifth BC bvp4c encounters a singular jacobian. Without the fifth BC or theta as a variable parameter, the only control over the launch angle would be from within the function, but I will be recursively changing theta to determine the angle with the lowest velocity after considering drag.
But besides from that, even without the fifth BC, or the variable parameter, the trejectory still launches at 89.96º every time and the whole point of including the variable parameter and extra BC is to try and control that.

Sign in to comment.

Accepted Answer

Max Friedman
Max Friedman on 1 Mar 2021
I ended up going a different route where I used ode45 to generate a hit point for a given initial velocity and angle, and a minimization function to change the velocity until it hit the correct target. This was nested inside of another minimization function which varried the angle to deduce the angle paired to the lowest velocity to hit the target. Once I had the optimal angle, it was simple enough to plug it back into the origional velocity calculation function to determine the required launch angle.
By using homemade bisection optimization functions, and fairly agressive velocity and angle errors of .1 degrees and m/s, it solves just about instantaneously.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!