hi there
i am trying to solve a simple problem using matlab. first, a bit of background for you:
a projectile is launched from a canon. after 15 seconds a parachute opens, the aim is to land the projectile exactly 10km away. various constants for drag, mass, parachute area etc are given in the problem. the aim of the excerise is to be able to 'put in' a desired distance at which to land the projectile and have matlab calculate the required angle to fire the canon at

I have written a code so far that uses the runge kutta method to solve the ODE, such that putting in the initial values and an angle will tell you where the projectile will land.

the code so far consists of 3 functions. it is run using the following command in matlab: Solve_Trajectory(0.1,0,0,900,45)  where the 45 is the angle of launch and can be varied. the code looks like this:


function 1:

Calculate the acceleration 

function [xdd, ydd] = find_acceleration(t, x, y, xd, yd, theta)



P = 1.207;                           % Air Density at sea level (kg/m^3)

if t <= 15, Apr = 0.25;   
else Apr = 1.1;                      % Projectile frontal area (m^2) (changes at t=15 when parachute opens)
end
   
M = 600;                             % Mass (kg)
   
g = 9.81;                            % Gravity

if t <=15, Cd = 0.38;                   
else Cd= 0.9;                        % Projectile drag coefficient (changes at t=15 when parachute opens)
end

Velocity  = ((xd^2)+(yd^2));         % Veloctiy squared (for drag)

Drag = 0.5*P*Velocity*Apr*Cd;        % Drag


% Calculations used to calculate xdd and ydd

xdd = (-Drag*cos(theta))/M;
ydd = ((-Drag*sin(theta))-M*g)/M;


function 2:
% Move the problem one step forward in time using Runge Kutta Method

function [x, xd, xdd, y, yd, ydd, theta] = Runge_kutta(t, x, xd, y, yd, theta, tInc)

%Calculate current acceleration

[xdd, ydd] = find_acceleration(t, x, y, xd, yd, theta);


%Define A,B,C,D, Beta and Delta for both x and y

%Ay and Ax

Ax=(tInc/2)*xdd;
Ay=(tInc/2)*ydd;

%Betay and Betax

Betax=(tInc/2)*(xd+(Ax/2));
Betay=(tInc/2)*(yd+(Ay/2));

[xdd, ydd] = find_acceleration(t+(tInc/2), x+Betax, y+Betay, xd+Ax, yd+Ay, theta);

%By and Bx

Bx=(tInc/2)*xdd;
By=(tInc/2)*ydd;

[xdd, ydd] = find_acceleration(t+(tInc/2), x+Betax, y+Betay, xd+Bx, yd+By, theta);

%Cy and Cx

Cx=(tInc/2)*xdd;
Cy=(tInc/2)*ydd;

%Deltay and Deltax

Deltax=tInc*(xd+Cx);
Deltay=tInc*(yd+Cy);

[xdd, ydd] = find_acceleration(t+(tInc), x+Deltax, y+Deltay, xd+2*Cx, yd+2*Cy, theta);

%Dy and Dx

Dx=(tInc/2)*xdd;
Dy=(tInc/2)*ydd;


%Calculate next set of values

x = x + tInc*(xd+(1/3)*(Ax+Bx+Cx));
xd = xd + (1/3)*(Ax+2*Bx+2*Cx+Dx);

y = y + tInc*(yd+(1/3)*(Ay+By+Cy));
yd = yd + (1/3)*(Ay+2*By+2*Cy+Dy);

theta = atan(yd/xd);

[xdd, ydd] = find_acceleration(t, x, y, xd, yd, theta);


function 3:

% Solve the initial-value ODE and plot the trajectory of the projectile

% inputs
% tInc = time increment (s)
% xAlpha = start horizontal displacement (m)
% yAlpha = start vertical displacement (m)
% Exit_V = exit velocity of the projectile at launch (m/s)
% angle = Angle of launch of the projectile (degrees)

function [t, x, y] = Solve_Trajectory(tInc, xAlpha, yAlpha, Exit_V, angle)

% Set initial conditions

theta(1) = angle*(pi/180);     % Convert launch angle into radians
t(1)=0;
x(1) = xAlpha;
xd(1) = Exit_V*cos(theta);
y(1) = yAlpha;
yd(1) = Exit_V*sin(theta);

% Loop through until projectile touches down
i=2;
while y >= 0
    [x(i), xd(i), xdd(i-1), y(i), yd(i), ydd(i-1), theta(i)] = Runge_kutta(t(i-1), x(i-1), xd(i-1), y(i-1), yd(i-1), theta(i-1), tInc);
    t(i)=t(i-1)+tInc;
    i=(i+1);
end

% Plot the results

plot(x, y);
title('Trajectory')
xlabel('Horizontal Distance [m]')
ylabel('Vertical Height [m]')

a = i-1;
x = x(a)


when running  Solve_Trajectory(0.1,0,0,900,45)  it runs through all of that and plots the trajectory of the projectile, and gives the distance travelled, x.

i appreciate thats a lot to look through, and if you're still reading then i thank you!

i now wish to write a shooting method solver to deduce what angle must be put in to the Solve_Trajectory function to output x=10000

can anyone help me with how to do this? my matlab skill are not great and it's taken me weeks to get this far!

many thanks if anyone can help
:)

> hi there
> i am trying to solve a simple problem using matlab.
> first, a bit of background for you:
> a projectile is launched from a canon. after 15
> seconds a parachute opens, the aim is to land the
> projectile exactly 10km away. various constants for
> drag, mass, parachute area etc are given in the
> problem. the aim of the excerise is to be able to
> 'put in' a desired distance at which to land the
> projectile and have matlab calculate the required
> angle to fire the canon at
> 
> I have written a code so far that uses the runge
> kutta method to solve the ODE, such that putting in
> the initial values and an angle will tell you where
> the projectile will land.
> 
> the code so far consists of 3 functions. it is run
> using the following command in matlab:
> Solve_Trajectory(0.1,0,0,900,45)  where the 45 is the
> angle of launch and can be varied. the code looks
> like this:
> 
> 
> function 1:
> 
> Calculate the acceleration 
> 
> function [xdd, ydd] = find_acceleration(t, x, y, xd,
> yd, theta)
> 
> 
> 
> P = 1.207;                           % Air Density at
> sea level (kg/m^3)
> 
> if t <= 15, Apr = 0.25;   
> else Apr = 1.1;                      % Projectile
> frontal area (m^2) (changes at t=15 when parachute
> opens)
> end
>    
> M = 600;                             % Mass (kg)
>    
> g = 9.81;                            % Gravity
> 
> if t <=15, Cd = 0.38;                   
> else Cd= 0.9;                        % Projectile
> drag coefficient (changes at t=15 when parachute
> opens)
> end
> 
> Velocity  = ((xd^2)+(yd^2));         % Veloctiy
> squared (for drag)
> 
> Drag = 0.5*P*Velocity*Apr*Cd;        % Drag
> 
> 
> % Calculations used to calculate xdd and ydd
> 
> xdd = (-Drag*cos(theta))/M;
> ydd = ((-Drag*sin(theta))-M*g)/M;
> 
> 
> function 2:
> % Move the problem one step forward in time using
> Runge Kutta Method
> 
> function [x, xd, xdd, y, yd, ydd, theta] =
> Runge_kutta(t, x, xd, y, yd, theta, tInc)
> 
> %Calculate current acceleration
> 
> [xdd, ydd] = find_acceleration(t, x, y, xd, yd,
> theta);
> 
> 
> %Define A,B,C,D, Beta and Delta for both x and y
> 
> %Ay and Ax
> 
> Ax=(tInc/2)*xdd;
> Ay=(tInc/2)*ydd;
> 
> %Betay and Betax
> 
> Betax=(tInc/2)*(xd+(Ax/2));
> Betay=(tInc/2)*(yd+(Ay/2));
> 
> [xdd, ydd] = find_acceleration(t+(tInc/2), x+Betax,
> y+Betay, xd+Ax, yd+Ay, theta);
> 
> %By and Bx
> 
> Bx=(tInc/2)*xdd;
> By=(tInc/2)*ydd;
> 
> [xdd, ydd] = find_acceleration(t+(tInc/2), x+Betax,
> y+Betay, xd+Bx, yd+By, theta);
> 
> %Cy and Cx
> 
> Cx=(tInc/2)*xdd;
> Cy=(tInc/2)*ydd;
> 
> %Deltay and Deltax
> 
> Deltax=tInc*(xd+Cx);
> Deltay=tInc*(yd+Cy);
> 
> [xdd, ydd] = find_acceleration(t+(tInc), x+Deltax,
> y+Deltay, xd+2*Cx, yd+2*Cy, theta);
> 
> %Dy and Dx
> 
> Dx=(tInc/2)*xdd;
> Dy=(tInc/2)*ydd;
> 
> 
> %Calculate next set of values
> 
> x = x + tInc*(xd+(1/3)*(Ax+Bx+Cx));
> xd = xd + (1/3)*(Ax+2*Bx+2*Cx+Dx);
> 
> y = y + tInc*(yd+(1/3)*(Ay+By+Cy));
> yd = yd + (1/3)*(Ay+2*By+2*Cy+Dy);
> 
> theta = atan(yd/xd);
> 
> [xdd, ydd] = find_acceleration(t, x, y, xd, yd,
> theta);
> 
> 
> function 3:
> 
> % Solve the initial-value ODE and plot the trajectory
> of the projectile
> 
> % inputs
> % tInc = time increment (s)
> % xAlpha = start horizontal displacement (m)
> % yAlpha = start vertical displacement (m)
> % Exit_V = exit velocity of the projectile at launch
> (m/s)
> % angle = Angle of launch of the projectile (degrees)
> 
> function [t, x, y] = Solve_Trajectory(tInc, xAlpha,
> yAlpha, Exit_V, angle)
> 
> % Set initial conditions
> 
> theta(1) = angle*(pi/180);     % Convert launch angle
> into radians
> t(1)=0;
> x(1) = xAlpha;
> xd(1) = Exit_V*cos(theta);
> y(1) = yAlpha;
> yd(1) = Exit_V*sin(theta);
> 
> % Loop through until projectile touches down
> i=2;
> while y >= 0
> [x(i), xd(i), xdd(i-1), y(i), yd(i), ydd(i-1),
> -1), theta(i)] = Runge_kutta(t(i-1), x(i-1), xd(i-1),
> y(i-1), yd(i-1), theta(i-1), tInc);
>     t(i)=t(i-1)+tInc;
>     i=(i+1);
> end
> 
> % Plot the results
> 
> plot(x, y);
> title('Trajectory')
> xlabel('Horizontal Distance [m]')
> ylabel('Vertical Height [m]')
> 
> a = i-1;
> x = x(a)
> 
> 
> when running  Solve_Trajectory(0.1,0,0,900,45)  it
> runs through all of that and plots the trajectory of
> the projectile, and gives the distance travelled, x.
> 
> i appreciate thats a lot to look through, and if
> you're still reading then i thank you!
> 
> i now wish to write a shooting method solver to
> deduce what angle must be put in to the
> Solve_Trajectory function to output x=10000
> 
> can anyone help me with how to do this? my matlab
> skill are not great and it's taken me weeks to get
> this far!
> 
> many thanks if anyone can help
> :)

Why don't you choose a ready-to-use MATLAB solver 
for your problem, e.g. BVP4C ?

Best wishes
Torsten.

hey
thanks for looking at this!

im not familiar with that im afraid, like i said my matlab skills are very poor!

anyway, i have to be able to demonstrate a shooting method solver for a given distance as part of the solution

thanks

"pat collett" wrote in message <icismi$2bj$1@fred.mathworks.com>...
> hey
> thanks for looking at this!
> 
> im not familiar with that im afraid, like i said my matlab skills are very poor!
> 
> anyway, i have to be able to demonstrate a shooting method solver for a given distance as part of the solution
> 
> thanks

hi, i know this tread is old, but how did you solve this problem?
thanks