Setting Boundary Conditions & Using bvp4c

Hi, I'm having trouble understanding how to use bvp4c. For my assignment, I have created a game that shoots a cannonball at a certain velocity and angle (which are the game inputs), it an attempt to hit a castle (the location is generated by the game).

Now, I need to make a cheat script using bvp4c that when given the position of the castle, outputs a set of velocity and angle values that would make the cannonball hit the castle. The castle is the red square in the image.

sol =  = bvp4c(@gamePhysics,@bcfun,solinit)

Here are the two second order ODE's that describe projectile's motion:

From there I converted them to these four first order ODE's:

function dydt = gamePhysics(t,y)
 m = 5;      % kg
 D = 0.013;      % N*s^2/m^2
 g = 9.81;      % m/s^2
dydt = zeros(4,1); 
dydt(1) = y(2)
dydt(2) = ((-D/m)*y(2)*(y(2)^2+y(4)^2)^(0.5));
dydt(3) = y(4);
dydt(4) = ((-g)-(D/m)*y(4)*(y(2)^2+y(4)^2)^(0.5));
end

I'm now stuck on making a function that describes the boundary conditions. At the origin, y(0) = 0 and at the castle, y(0) = x_castle.

Also, I'm confused on what the solinit function in bvp4c means. Am I guessing a random velocity and angle?

Thank you! Sorry, I'm pretty lost on how to use this function.

Answer 1

[Torsten]

xcastle = 20;                               % x-coordinate of target
v0 = 20;                                    % Prescribed (absolute) velocity at (0,0)
dxdt0 = 1;
dydt0 = sqrt(v0^2-dxdt0^2);                 % Initial guesses for x- and y-component of velocity at (0,0)
sol = fminunc(@(u)opt(u,xcastle,v0),dxdt0); % Call optimizer

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

dxdt = sol(1)

dxdt = 5.5300

dydt = sqrt(v0^2-dxdt^2)                    % Compute x- and y-component of velocity at (0,0)

dydt = 19.2203

theta = atan(dydt/dxdt)*180/pi              % Compute angle at (0,0)

theta = 73.9485

tspan = 0:0.05:1000;                        % plot trajectory
[~,t,y] = solve_ode(dxdt,dydt,tspan);
plot(y(:,1),y(:,3))
axis equal

% Minimize distance between x-coordinate of target and intersection of
% trajectory with ground level
function res = opt(u,xcastle,v0)
  dxdt = u(1);
  dydt = sqrt(v0^2-dxdt^2);
  tspan = [0 1000];
  xcut = solve_ode(dxdt,dydt,tspan);
  res = (xcut - xcastle)^2;
end

% Solve ode that defines the trajectory
% Stop integration when ground level is hit
function [xcut,t,y] = solve_ode(dxdt,dydt,tspan)
  u0 = [0 dxdt 0 dydt];
  options = odeset('Events',@Event);
  [t,y] = ode45(@gamePhysics,tspan,u0,options);
  xcut = y(end,1);
end

% Define the physical equations
function dydt = gamePhysics(t,y)
  m = 5;        % kg
  D = 0.013;    % kg/m  
  g = 9.81;     % m/s^2 
  dydt = zeros(4,1); 
  dydt(1) = y(2);
  dydt(2) = ((-D/m)*y(2)*(y(2)^2+y(4)^2)^(0.5));
  dydt(3) = y(4);
  dydt(4) = ((-g)-(D/m)*y(4)*(y(2)^2+y(4)^2)^(0.5));
end

% Stop integration when trajectory intersects the ground level
function [value,isterminal,direction] = Event(t,y) 
  value = y(3);
  isterminal = 1;
  direction = -1;
end