Shooting method to solve ODE problem 1

[boileroo]

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

 

[GregLocock]

Given the actual non linearity in the system it seems unlikely that there is an analytical solution. Therefore it seems likely that an intelligent guessing algorithm is required. Mr Newton and Mr Raphaelson might have written one.

Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[FeX32]

yea, ODE45 is very useful for that


Fe

 

[electricpete]

If I read the op right, a simulation has already been implemented equivalent to ode45

In that case, all that is required is to build a loop to do multiple simulations varying theta each time and plot final position vs theta. Look for value of theta where the curve crosses the target position.

boileroo - you must already know how to build a loop and generate a plot because you did it in your posted code. Can't you just build one more loop and generate one more plot as I described above? What part of that process are you having trouble with.

=====================================
(2B)+(2B)' ?

 

[SomptingGuy]

The "Shooting Method" is well defined:

http://en.wikipedia.org/wiki/Shooting_method

Numerical recipes also describes it very well.

It doesn't involve time-stepping the ODEs.

- Steve

 

[electricpete]

I'm not sure what you are talking about. The shooting method is used when the boundary conditions define the problem but do not provide a full set of initial conditions for initial value problems. The unknown initial conditions are varied to match the known non-initial-condition boudary value. It is well suited for numeric solution on problems such as this. One of us is mistaken I think.... if it is me, please explain.

=====================================
(2B)+(2B)' ?

 

[SomptingGuy]

I need to re-read my NRC. Best book ever for describing algorithms and their benefits/limitations in human language. Not very good for real code.

- Steve

 

[hacksaw]

Simple, for launch angles >45 Deg, pick a launch angle, shoot, if past the target raise the launch angle. Below 45 Deg the logic reverses. Propellant charge adjustment follows similar logic.

 

[FeX32]

caboom


Fe