iconEuler Examples

Ballistic Shots

by R. Grothmann

We want to demonstrate the use of numerical analysis with Euler for
real world problems.

The problem is to shoot a canon ball with speed v and angle a. This
can, of course, be computed exactly. But first we restrict ourselves
to numerical solutions.

We fix the speed v=10 m/s and approximate the gravition with g=10
m/s^2.
>v=10; g=10;
Now we plot a specific ballistic curve for a canon ball fired at 45°
with the well known formulas for the curve neglecting friction.

Time runs from 0 s to 2 s.
>t=0:0.1:2;
The movement in x-direction (height) is cos(a)*v*t. The movement in
y-direction is

Ballistic Shots

>plot2d(cos(45°)*v*t,sin(45°)*v*t-g*t^2/2):

Ballistic Shots

We notice that we are not interested in the negative y-values and clip
the plot to the positive quadrant.
>plot2d(cos(45°)*v*t,sin(45°)*v*t-g*t^2/2,a=0,b=10,c=0,d=10):

Ballistic Shots

We now compute, when the ball has its highest altitude. We use the
function fmax, which computes the maximum of another function or
expression in a given range.
>tmax=fmax("sin(45°)*v*x-g*x^2/2",0,2)
0.707106772995
How high is the maximal height?
>x=tmax; sin(45°)*v*x-g*x^2/2
2.5
To make things a little bit more comfortable, we introduce two
functions computing the point (sx,sy) at time t.

We take v and g from the global values, but add a as a parameter to
these functions.
>function sx (t,a) := cos(a)*v*t;
>function sy (t,a) := sin(a)*v*t-g*t^2/2;
This makes computing the maximal height easier.
>a=45°; sy(fmax("sy(x,a)",0,2),a)
2.5
We wish to compare several angles. Thus we generate a vector of angles
from 5° to 85° and generate all curves for all these angles
simultanously.

To do so, we need the angle a as a column vector and the time t as a
row vector. Combining both yields a matrix, with a curve in each row
for each angle.
>a=(5:5:85)°; a=a';
>plot2d(sx(t,a),sy(t,a),a=0,b=10,c=0,d=10):

Ballistic Shots

Lets us now compute the distance from the origin to the contact of the
ball with the ground. We use the simplest of methods, the bijection
method. Since 0 is already a solution of sy=0, we start a little bit
off.

We seek the zero of vy. The additional parameter a (after the
semicolon) is passed from bisect to sy.

Note: Our function "impact" will only work for scalar a, not for
vectors a, if we do not map.
>function map distance (a) := sx(bisect("sy",0.00001,2;a),a)
What happens, if we take 50°?
>distance(50°)
9.84807753012
We get about 9.85m.

Let us plot the impact distance with varying angles.
>a=(1:1:89)°;
To compute the impact function for all angles, we apply "distance" to
the vector a. The function is mapped to all elements of a.
>plot2d(deg(a),map("distance",a)):

Ballistic Shots

We can also compute the optimal angle and its distance.
>amax=fmax("distance",1°,89°); deg(amax), distance(amax),
45.0000003104
10
Of course, it is easy to invert the function sx, i.e., compute
the time t for the value x. We do not need to do that numerically.

So we write a function that computes the height of the ball
at distance x, if the angle is a.
>function height (x,a) := sy(x/(cos(a)*v),a)
Test this with 45° from 0 to 10.
>plot2d("height";45°,a=0,b=10);
We now want to get as high as possible in 5m distance. Again,
we use fmax, but the argument x is now the second parameter to
the function height, the angle.
>degprint(fmax("height(5,x)",1°,89°))
63°26'5.82''
We need the angle argument as first argument for the next function.
>function height1(a,x) := height(x,a)
Now, we make a function that computes the angles that yields the
maximal heights for given distances. Since fmax calls height to
maximize the angle, we call height1(x,xx).
>function heighestangle (xx) := fmax("height1",1°,89°;xx)
Compute various maximal heights and plot them together with
some shots.
>x=0.2:0.05:10; y=map("heighestangle",x);
>plot2d(x,height(x,y),color=red,thickness=2);
>a=(5°:5°:85°)'; plot2d(sx(t,a),sy(t,a),add=1):

Ballistic Shots

Is the curve of maximal heights a parabola?

Yes it is. To demonstrate, we compute the interpolation polynomial of
second degree through the know values at -10, 0 and 10.
>d=interp([-10,0,10],[0,5,0]);
Then we plot the polynomial over the current view in color 13.
Indeed, we get the same curve.
>plot2d("interpval([-10,0,10],d,x)",add=1,color=10,thickness=3):

Ballistic Shots

Next, we ask how to reach the point x=2, y=2. There seem to
be two solutions. We get both by searching the angle below and
the angle above the one with maximal height.
>a=heighestangle(2); degprint(a)
78°41'24.24''
>height(2,a)
4.8
>a1=bisect("height(2,x)-2",0,a); degprint(a1)
51°31'33.49''
>a2=bisect("height(2,x)-2",a,pi/2); degprint(a2)
83°28'26.51''
Let's plot both solutions to one plot.
>t=0:0.01:2;
>plot2d(sx(t,a1),sy(t,a1),a=0,b=5,c=0,d=5);
>plot2d(sx(t,a2),sy(t,a2),add=1):

Ballistic Shots

Assume we shoot a ball in 80° and another one in 40°. Can we hit the
first with the second in the air?

We plot the two paths.
>plot2d("height(x,80°)",a=0,b=4); plot2d("height(x,40°)",add=1):

Ballistic Shots

Now we compute the point, where both paths meet and compute
the time difference.
>xd=bisect("height(x,80°)-height(x,40°)",0.001,5)
3.07201661702
>xd/(v*cos(80°))-xd/(v*cos(40°))
1.3680805733

Friction

Now we introduce friction. We assume the friction force is
proportional to c*v^2, where v is the speed.

We need to solve a differential equation, using the "runge" function.
The state x of the ball at any time t is 4-dimensional and contains
two point coordinates and two speed coordinates. We need a function to
compute x'=f(t,x).

We set x=[sx,sy,sx',sy'].
>function f(t,x,c=0.005) ...
global g;
v=sqrt(x[3]*x[3]+x[4]*x[4]); rho=c*v;
return [x[3],x[4],-rho*x[3],-g-rho*x[4]]
endfunction
We need to pass a vector of times t, where x will be computed.
"runge" returns a matrix s, with one value for s in each column.
We may plot the path of the ball from the first two rows of
s.
>t=0:0.01:2; s=ode("f",t,[0,0,cos(45°)*v,sin(45°)*v]);
>plot2d(s[1],s[2],a=0,b=10,c=0,d=10);
As you see, we get a little shorter.
>plot2d(sx(t,45°),sy(t,45°),add=1,color=blue):

Ballistic Shots

Now we look at the speed over a larger time. The speed will
tend to some value, where friction is equal to gravity.
>t=0:0.01:20; s=ode("f",t,[0,0,cos(45°)*v,sin(45°)*v]);
>plot2d(t,sqrt(s[3]^2+s[4]^2)):

Ballistic Shots

We wish to compute the angle of maximal width for this problem.
It will no longer be exactly 45°.

First we define a function returning the height at time t.
>function height (t,a) ...
global v;
s=runge("f",[0,t],[0,0,cos(a)*v,sin(a)*v],20);
return s[2,2];
endfunction
>plot2d("height";45°,a=0,b=2):

Ballistic Shots

Then we seek the time when the height is 0 for any angle a,
and return the width at that time.
>function width (a) ...
t=bisect("height",0.001,2;a);
global v;
s=runge("f",[0,t],[0,0,cos(a)*v,sin(a)*v],20);
return s[1,2];
endfunction
>width(45°)
9.62578470286
Now we maximize the width over the angles.

This computation takes some time (approx. 3 seconds on my modern
Notebook)

We find that the optimal angle is slightly less than 45°.
>amax=fmax("width",30°,60°); degprint(amax)
44°42'25.31''

Symbolic Solution

Now we compute ballistic shots with Maxima, with and without air
resistance.

First the x-coordinate of the shot path (neglecting friction).
>function sx(v,t,a) &= v*t*cos(a)
                              cos(a) t v

For the y-coordinate we subtract the influence of gravity.
>function sy(v,t,a) &= v*t*sin(a)-g*t^2/2
                                          2
                                       g t
                          sin(a) t v - ----
                                        2

Let us plot a sample shot.
>v=10; g=10; a=%pi/4; plot2d("sx(v,x,a)","sy(v,x,a)", ...
 a=0,b=10,c=0,d=10,xmin=0,xmax=2):

Ballistic Shots

Now we solve for the time, when the ball hits the ground.
>&solve(sy(v,t,a)=0,t), function wa(v,a) &= at(sx(v,t,a),%[1])
                            2 sin(a) v
                       [t = ----------, t = 0]
                                g


                                           2
                          2 cos(a) sin(a) v
                          ------------------
                                  g

Let us plot the width with varying angles.
>plot2d("wa(v,rad(x))",0,90):

Ballistic Shots

We can compute the angle of maximum distance numerically in Euler.
>deg(fmax("wa(v,x)",0,90°))
44.9999996606
The same in Maxima using calculus.
>&solve(diff(wa(v,a),a)=0,a)
                 [sin(a) = - cos(a), sin(a) = cos(a)]

We need to help Maxima.
>&solve(trigreduce(diff(wa(v,a),a))=0,a)
                                    pi
                               [a = --]
                                    4

Let us compute the height of a ball at distance x.
>&solve(sx(v,t,a)=x,t), function hx(x,v,a) &= ratsimp(at(sy(v,t,a),%))
                                    x
                            [t = --------]
                                 cos(a) v


                                      2        2
                     2 cos(a) sin(a) v  x - g x
                     ---------------------------
                                 2     2
                            2 cos (a) v

Plot a shot at angle 45°.
>plot2d("hx(x,v,45°)",a=0,b=10,c=0,d=10):

Ballistic Shots

Now we use the Euler matrix language to plot shots for angles
from 45° to 90°.
>x:=0:0.01:10; y:=hx(x,v,(45°:5°:90°)'); ...
 plot2d(x,y,a=0,b=10,c=0,d=10):

Ballistic Shots

We wish to compute the resultant. In this case, it the maximal
height one can reach at distance x.
>&trigsimp(diff(hx(x,v,a),a)), &solve(%=0,a)
                              2               2
                      cos(a) v  x - sin(a) g x
                      -------------------------
                                3     2
                             cos (a) v


                                           2
                                   cos(a) v
                         [sin(a) = ---------]
                                      g x

This is not quite complete, but we can use it nevertheless.
>function par(x,v) &= ratexpand(at(hx(x,v,a),a=expand(atan(v^2/(g*x)))))
                               2       2
                              v     g x
                              --- - ----
                              2 g      2
                                    2 v

We add this parabola to the plot above.
>plot2d("par(x,v)",add=1,color=13,thickness=3):

Ballistic Shots

This does also answer the question, how far we can get if we aim
from some level b different from 0.

E.g., assume we are b meters above ground and throw a stone with
speed v at an optimal angle, at gravity g. The distance is then
our resolvant at y=-b.
>&solve(par(x,v)=-b,x)[2]
                                    2
                            v sqrt(v  + 2 b g)
                        x = ------------------
                                    g

To find the optimal angle for such a throw, we need tu use numerical
methods. First we find the distance, we can reach, if we use an angle
a.
>&solve(sy(v,t,a)=-b,t), function wa(a,v) &= at(sx(v,t,a),%[2])
                                2     2
             sin(a) v - sqrt(sin (a) v  + 2 b g)
        [t = -----------------------------------, 
                              g
                                          2     2
                                  sqrt(sin (a) v  + 2 b g) + sin(a) v
                              t = -----------------------------------]
                                                   g


                              2     2
            cos(a) v (sqrt(sin (a) v  + 2 b g) + sin(a) v)
            ----------------------------------------------
                                  g

We use the Euler fmax function to find the angle for the maximal
distance at height 5, with v=10, g=10 as above.
>v:=10; g:=10; b:=5; deg(fmax("wa(x,v)",0,90°))
35.2643898543

Friction

Let us try to add some air resistance. However, we can only solve that
equation exactly, if we assume the air resistance is quadratic with
the speed and that we simply drop the body.

We set up a differential equation for the speed at time x.
>eq1 &= "'diff(v(x),x) = -g+k*v(x)^2"
                       d              2
                       -- (v(x)) = k v (x) - g
                       dx

This can be solved.
>&assume(g>0); &assume(k>0); sol &= ode2(eq1,v(x),x)
                    k v(x) - sqrt(g) sqrt(k)
                log(------------------------)
                    k v(x) + sqrt(g) sqrt(k)
                ----------------------------- = x + %c
                      2 sqrt(g) sqrt(k)

The solution is implicit. We need to solve for v(x) first, and
assign the solution to sol(x).
>function sol(x) &= ratsimp(rhs(solve(sol,v(x))[1]))
        
                       2 sqrt(g) sqrt(k) x + 2 %c sqrt(g) sqrt(k)
     sqrt(g) sqrt(k) (E                                           + 1)
   - -----------------------------------------------------------------
                2 sqrt(g) sqrt(k) x + 2 %c sqrt(g) sqrt(k)
             k E                                           - k

Next we determine the constant %c, so that the speed at x=0
is 0.
>constants &= solve(sol(0)=0,%c)
                      log(- I)               log(I)
             [%c = ---------------, %c = ---------------]
                   sqrt(g) sqrt(k)       sqrt(g) sqrt(k)

Substitute the solution into sol(x) and fix sol(x).
>function solution(x) &= at(sol(x),constants[2])
                                     2 sqrt(g) sqrt(k) x
               sqrt(g) sqrt(k) (1 - E                   )
             - ------------------------------------------
                           2 sqrt(g) sqrt(k) x
                      - k E                    - k

Next we plot the speed versus the time. As expected, we get
a limit speed.
>k:=0.01; plot2d("solution(x)",0,40):

Ballistic Shots

We can find this limit exactly with Maxima.
>&limit(solution(x),x,inf)|k=1/100|g=10, %()
                                   3/2
                               - 10

-31.6227766017
Now we integrate the speed and find the drop length at time
a.
>&assume(x>0); function h(x) &= integrate(solution(t),t,0,x)
                2 sqrt(g) sqrt(k) x
           log(E                    + 1)   sqrt(g) x   log(2)
         - ----------------------------- + --------- + ------
                         k                  sqrt(k)      k

If we plot this, it looks like a parabola.
>plot2d("h(x)",0,5):

Ballistic Shots

With k=0, this converges to a parabola. We need to compute the limit
using Taylor a series expansion.
>&tlimit(h(x),k,0)
                                     2
                                  g x
                                - ----
                                   2

Examples Homepage