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

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

y=deval(sol,linspace(1,d,100));
plot(y(3,:),y(1,:));
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);
end

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));
end

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;
end

[darova]

Why acceleration depends on ?

[Max Friedman]

those are the drag components, and drag is dependent on v^2. k is defined earlier as 1/2 * rho * C_d * A / m. Though this shouldn't come into play as C_d is 0.

[darova]

You have too much boundary conditions. There is no need of angle since you have vectors (vx,vy)

4 equations - 4 conditoins

[Max Friedman]

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.

Answer 1

[Max Friedman]

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.