For a uni project I was asked to code a program to solve a two-point boundary value problem. The problem is referred as the putting problem and illustrates the shooting method for two-point boundary-value problems. The differential equations are developed and numerical results are given.
My differential system is: y3*y5, y4*y5, g*n3*(n1 - mu*y3/s)*y5, g*n3*(n2 - mu*y4/s)*y5, 0
To solve the differential equations developed using newtons laws I coded a 4th order Runge-Kutta method for a differential system. My problem is that now I have to approximante the initial speed that has to be given to the ball and also the time T that it took the ball to reach the hole (or y3(0), y4(0) and y5(0) that I will respectively name a,b and T). For that I find useful to create a residual system that needs to be equalized to 0 using Newton-Raphson. To clarify, at t = 1 (or when the ball is supposed to fall in the hole) y1(1) - 20 needs to be 0, y2(2) - 0 also 0 and y3^2(1) + y4^2(1) + x'3^2(1) = 0 (the first 2 being the position and the 3rd the velocity) assuming that the position of the hole is (20, 0). To apply Newton-Raphson I need a equation system however, here, I have only numerical results after applying RK4. Basicaly my code has to iterate between RK4 and Newton-Raphson with each time, better values of a,b and T. Of course the residual system can't be exactly 0 so we can take an error of 0.0015 for both the speed and the position.
This is my code with the boundary conditions of a = 4, b = 0 and T = 5 with a solution using Runge-Kutta:
import matplotlib.pyplot as plt
import numpy as np
#import numpy as np
from sympy import Symbol, Function, diff
def surface_plane(x1, x2):
return 0
def surface_parabolique(x1, x2):
return ((x1-10)**2)/125 + ((x2-5)**2)/125-1
def surface_inclinee(x1,x2, alpha=0):
return alpha*x2
def F(t_, Y, mu=.2,g = 9.81):
""" Attention, la dependance en t_ est implicite."""
# surface est une fonction globale
# les differentes variables
y1,y2,y3,y4,y5 = Y
## evaluation de la normale a la surface
normale = np.array(( diff(surface,x1).subs([(x1,y1),(x2,y2)]),
diff(surface,x2).subs([(x1,y1),(x2,y2)]),
1)).astype(float)
#print("normale :",normale,end="\t")
n1 = -normale[0]/np.linalg.norm(normale)
n2 = -normale[1]/np.linalg.norm(normale)
n3 = normale[2]/np.linalg.norm(normale)
# vitesse : evaluation de la composante verticale
# -----------------------------------------------
if isinstance(v3,(float,int)):
v3_ = v3
else:
# evaluation de la derivee
v3_ = v3.subs([(x1_,y1),(diff(x1_),y3),
(x2_,y2),(diff(x2_),y4)]).doit()
#print("v3_ :",v3_)
# vitesse : evaluation de la norme
# --------------------------------
global s
s = (y3**2+y4**2+v3_**2)**.5
#fG =
return np.array([y3*y5,
y4*y5,
g*n3*(n1 - mu*y3/s)*y5,
g*n3*(n2 - mu*y4/s)*y5,
0])
def rk4(F, t0, y0, h, n):
y = np.zeros((n+1,len(y0)))
y[0,:] = y0
for i in range(n):
k1 = h * F(t0 +i*h , y[i,:] )
k2 = h * F(t0 +(i+.5)*h, y[i,:] + 0.5*k1)
k3 = h * F(t0 +(i+.5)*h, y[i,:] + 0.5*k2)
k4 = h * F(t0 +(i+1)*h , y[i,:] + k3 )
y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
return y
if __name__ == "__main__":
mu = 0.2
t = Symbol('t')
x1 = Symbol('x1')
x1_ = Function('x1')(t)
x2 = Symbol('x2')
x2_ = Function('x2')(t)
#surface = surface_plane(x1,x2)
surface = surface_parabolique(x1,x2)
#surface = surface_inclinee(x1,x2)
print("surface :",surface)
if isinstance(surface,(float,int)):
v3 = surface
# on a une constante
else:
#print("surface.subs([(x1,x1_),(x2,x2_)])",surface.subs([(x1,x1_),(x2,x2_)]))
# on a une fonction de t que l'on derive
v3 = diff(surface.subs([(x1,x1_), (x2,x2_)]), t)
print("v3 :",v3)
t0 = 0 # temps initial
n = 100 # pas pour RK4
tau = np.linspace(t0, 1, n+1)
h = tau[1]-tau[0] # dt ou h pour aller de 0 a 1
a = 4 # v1 initiale
b = 0 # v2 initiale
T = 5 # tps d'arrivee
# conditions initiales
Y0=(0, 0, a, b, T)
# on applique runge-kutta d'ordre 4
yn = rk4(F, t0, Y0, h, n)
print(yn[(0,-1),:])
plt.plot(tau,yn,'.')
plt.legend(('y1','y2','y3','y4','y5'))
