laddersSoluce.py
# -------------------------------------------------------------------------
#
# PYTHON for DUMMIES
# Probleme 9
#
# Solution detaillée
# Vincent Legat
# Nathan Coppin
#
# -------------------------------------------------------------------------
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
# -------------------------------------------------------------------------