shoot.py
# -------------------------------------------------------------------------
#
# PYTHON for DUMMIES 19-20
#
# Shooting technique
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#
# -------------------------------------------------------------------------
from numpy import *
from matplotlib import pyplot as plt
import matplotlib
matplotlib.rcParams['toolbar'] = 'None'
# plt.rcParams['figure.facecolor'] = 'lavender'
# plt.rcParams['axes.facecolor'] = 'lavender'
# -------------------------------------------------------------------------
#
# -1- Intégration de la trajectoire de l'obus
#
# dxdt(t) = u(t)
# dudt(t) = - gamma*u(t)
# dydt(t) = v(t)
# dvdt(t) = - 9.81 - gamma*v(t)
#
#
def f(u):
friction = 0.1 * sqrt(u[1]*u[1]+u[3]*u[3])
dxdt = u[1]
dudt = - friction *u[1]
dydt = u[3]
dvdt = -9.81 - friction * u[3]
return array([dxdt,dudt,dydt,dvdt])
# -------------------------------------------------------------------------
#
# -2- Schema de Runge-Kutta classique d'ordre 4
#
def rungekutta(a,h):
imax = int(8/h)
u = a * cos(pi/4)
v = a * sin(pi/4)
X = arange(imax+1)*h
U = zeros((imax+1,4)); U[0,:] = [0,u,0,v]
i = 0
while (U[i,2] >= 0) and (i < imax) :
K1 = f(U[i,:] )
K2 = f(U[i,:]+K1*h/2)
K3 = f(U[i,:]+K2*h/2)
K4 = f(U[i,:]+K3*h )
U[i+1,:] = U[i,:] + h*(K1+2*K2+2*K3+K4)/6
i = i+1
return i,X,U
#
# -3- Méthode de bissection
#
def shoot(a,h):
i,X,U = rungekutta(a,h)
plt.plot(U[:i,0],U[:i,2],'-b')
return (35-U[i,0])
def shootingMethod(h,tol):
n = 1; nmax = 40;
a = 100; fa = shoot(a,h)
b = 2000; fb = shoot(b,h)
n = 0; delta = (b-a)/2; x = a;
if (fa*fb > 0) :
raise RuntimeError('Bad initial interval')
while (abs(delta) >= tol and n <= nmax) :
delta = (b-a)/2; n = n + 1;
x = a + delta; fx = shoot(x,h)
print(" x = %14.7e (Estimated error %13.7e at iteration %d)" % (x,abs(delta),n))
if (fx*fa > 0) :
a = x; fa = fx
else :
b = x; fb = fx
if (n > nmax) :
raise RuntimeError('Too much iterations')
return x
# -------------------------------------------------------------------------
fig = plt.figure("Shooting techniques")
h = 0.01
tol = 1e-1
a = shootingMethod(h,tol)
i,X,U = rungekutta(a,h)
plt.plot(U[:i,0],U[:i,2],'-r')
plt.show()
print(" ============ Final value for x(%f) = %.4f " % (X[i],U[i,0]))