×
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Are you an
Engineering professional?
Join Eng-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Shooting method to solve ODE problem

Shooting method to solve ODE problem

(OP)
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
:)
Replies continue below

RE: Shooting method to solve ODE problem

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: Eng-Tips.com Forum Policies  http://eng-tips.com/market.cfm?

RE: Shooting method to solve ODE problem

yea, ODE45 is very useful for that

Fe

RE: Shooting method to solve ODE problem

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)'  ?

RE: Shooting method to solve ODE problem

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)'  ?

RE: Shooting method to solve ODE problem

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

RE: Shooting method to solve ODE problem

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.

RE: Shooting method to solve ODE problem

caboom

Fe

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Close Box

Join Eng-Tips® Today!

Join your peers on the Internet's largest technical engineering professional community.
It's easy to join and it's free.

Here's Why Members Love Eng-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!