laddersSoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for DUMMIES 19-20
# Probleme 6
#
# Solution detaillée
#  Vincent Legat
#
# -------------------------------------------------------------------------




from scipy.linalg import norm,solve


# -------------------------------------------------------------------------
#
# -1- Schéma de Newton-Raphson pour les équations suivantes
#     obtenues avec des triangles semblables en élevant au carré les rapports !
#
#
#       a^2 x^2 = (x^2+c^2) (x^2 + y^2)
#       b^2 b^2 = (x^2+c^2) (x^2 + y^2)    
# 
#

def f(u):
  global a,b,c
  x = u[0]; y = u[1]
  return [a*a*x*x -(x*x+c*c)*((x+y)**2),
          b*b*y*y -(y*y+c*c)*((x+y)**2)]

def dfdx(u):
  global a,b,c
  x = u[0]; y = u[1]
  return [[2*a*a*x - 2*x*((x+y)**2) - 2*(x*x+c*c)*(x+y) ,       -2*(x*x+c*c)*(x+y) ],
          [-2*(y*y+c*c)*(x+y) ,      2*b*b*y - 2*y*((x+y)**2) - 2*(y*y+c*c)*(x+y) ]]


# -------------------------------------------------------------------------

def laddersSolve(geometry,tol,nmax):
  global a,b,c
  a = geometry[0]
  b = geometry[1]
  c = geometry[2]

#
# Astuce éventuelle pour le cas a=b qui rend la jacobienne de Newton-Raphson singulière
# Dans ce cas précis, on peut obtenir directement une solution analytique
# 
# On utilise la tolérance pour le test du cas particulier
# Dans le cas d'une itération simple, il a été nécessaire 
# d'intéger une tolérance absolue (beurck :-)
#
# 
#   if (abs(a-b) < tol) :
#     x = ((b*b - 4*c*c)**(1/2))/2
#     return solve([[1,0],[0,1]],[x,x])
# 
# En pratique, solve fournira la bonne réponse même sans y ajouter l'astuce
#

#
# Heuristique de Nathan pour obtenir un bon candidat initial
# Cela fonctionne remarquablement bien :-)
#

  r = (min(a,b) - 0.1)/2
  x = [r,r]  
  
  n = 0; delta = tol+1
  while (norm(delta) > tol and n < nmax):
    delta = - solve(dfdx(x),f(x))
    x = x + delta
    n = n+1
    print(" Estimated error %9.2e at iteration %d : " % (norm(delta),n),x)
  return x
  
# -------------------------------------------------------------------------

def laddersIterate(geometry,x):
  global a,b,c
  a = geometry[0]
  b = geometry[1]
  c = geometry[2]

  
  delta = - solve(dfdx(x),f(x))
  x = x + delta
  return x
  
  # -------------------------------------------------------------------------