I felt like sharing the solution of a simple yet interesting problem from the end-sem exam of our Numerical Analysis course. I’ve used R (learnt another language, yay!) and OriginPro. The solution involves solving differential equations and interpolation.
Edit: The position values in x,y,z have been given in kilometres.
Let’s break down the problem:
1. We have been given velocity equations of the plane, hence we can integrate them using the given initial position (0, 0, 0) and get the trajectory of the plane till the end of time. We can determine exactly where the plane is after 5 minutes in x, y, z coordinates. This is the time when it crosses over to the enemy country and becomes visible to the enemy radar.
2. Now we switch over to the enemy observer. He can observe the plane in polar coordinates (r, θ1, θ2). We only need to worry about r, as the missile launcher can calculate angles by itself. So we begin taking down observations for r each second. Although the value of r is readily available to the observer, we need to calculate it using the instantaneous position of the plane and position of the observer (simple distance between two points).
3. We wait till the plane enters into the sphere of 500m radius around the observer. Now the observer has to calculate where the plane will be 1 second later, and fire his missile accordingly with suitable speed. (Note that we know exactly where the plane is going, but the observer doesn’t. We cannot use our knowledge of the plane’s equation of motion.)
4. The observer needs to use the data he recorded in the past few seconds and predict a point where the plane is expected to be 1 second later.
Solution:
We use Runge-Kutta algorithm to integrate the velocity equations and obtain the values of position at each second. We start tracking the distance between plane and observer after 5 mins, till the distance becomes 500m.
[sourcecode title=”Track.R” wraplines=”false”]
library("pracma")
f <- function(t, s)
{
s1 <- 0.2*exp(-0.01*(t-220))/(1+exp(-0.01*(t-220)))^2
s2 <- 0.01994711*exp(-0.5*((t-200)/100)^2)
s3 <- 0.0052*exp(-t/100)
return(c(s1, s2, s3))
}
s0 <- c(0,0,0)
a <- 0; b <- 500
n <- b
sol <- rk4sys(f, a, b, s0, n)
time <- 301
radar <- matrix(nrow=10,ncol=2)
r <- 1000
while(r>500)
{
r = sqrt((sol$y[time,1]-12)^2+(sol$y[time,2]-4.15)^2+(sol$y[time,3]^2))
r = r*1000
time = time + 1
radar[time-301,1] <- time-302
radar[time-301,2] <- r
}
[/sourcecode]
The following file contains the calculated values of position (x, y, z) at each second starting from 0 seconds: Planepos.txt
Calculated values of r (time starting after 5 minutes):
0 s ——– 534.4396 m
1 s ——– 519.4016 m
2 s ——– 507.8468 m
3 s ——– 499.9667 m
Now we need to find a polynomial function of order 3 which fits these 4 points. It can be done using various algorithms for numerical interpolation. We then use the polynomial to find a 5th point, which will be the approximate position of the plane after 1 second (i.e. the 4th second).
Here I’ve used OriginPro to directly find a fitting polynomial after plotting the above points. However, this is just an approximate fit. The accuracy can be greatly improved with a code for interpolation.
The equation turns out to be r = 539.784 – 26.51383 t + 6.99025 t2 – 0.85882 t3.
Using the above equation, the calculated value of r at the 4th second = 490.6082 m.
Therefore, the missile needs to be fired after tracking the plane for 3 s, with a speed of 490.6082 m/s.