Ejercicio
Referencia: Chapra 6.5 p162, Chapra Ejercicio 6.11 p166
Con el método de Newton-Raphson para múltiples ecuaciones, determine las raíces para:
x^2+xy =10 y + 3xy^2 = 57Observe que un par correcto de raíces es x=2 y y=3.
Use como valores iniciales x=1.5, y=3.5
Planteamiento
Las ecuaciones se expresan de la forma f(x,y) = 0
x^2+xy -10 = 0 y + 3xy^2 -57 = 0Se puede usar extensiones de los métodos abiertos para resolver ecuaciones simples, por ejemplo Newton-Raphson.
u_{i+1} = u_i + (x_{i+1}-x_i)\frac{\partial u_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial u_i}{\partial y} v_{i+1} = v_i + (x_{i+1}-x_i)\frac{\partial v_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial v_i}{\partial y}ecuaciones que se pueden reordenar y encontrar la solución a partir de la matriz Jacobiano.
Algoritmo en Python
Usando un algoritmo para resolver el Jacobiano y estimar los puntos luego de cada iteración se obtienen:
iteración: 1
Jacobiano con puntos iniciales:
[ 6.5 1.5 ]
[ ]
[36.75 32.5]
determinante: 156.12499999999994
puntos xi,yi: 2.03602882305845 2.84387510008006
error: 0.656124899919936
iteración: 2
Jacobiano con puntos iniciales:
[6.91593274619696 2.03602882305845]
[ ]
[24.2628767545662 35.7412700376474]
determinante: 197.78430344142245
puntos xi,yi: 1.99870060905582 3.00228856292451
error: 0.158413462844444
iteración: 3
Jacobiano con puntos iniciales:
[6.99968978103616 1.99870060905582]
[ ]
[27.0412098452019 37.0040558756713]
determinante: 204.96962918261596
puntos xi,yi: 1.99999998387626 2.99999941338891
error: 0.00228914953559523
iteración: 4
Jacobiano con puntos iniciales:
[6.99999938114143 1.99999998387626]
[ ]
[26.9999894410015 36.9999926704397]
determinante: 204.9999473486533
puntos xi,yi: 1.99999999999998 3.00000000000008
error: 5.86611161867978e-7
Resultado:
1.99999999999998 3.00000000000008
Algoritmo presentado para dos ecuaciones y dos incógnitas, en la unidad 3 se puede ampliar la propuesta. Revisar el método de Gauss-Seidel y Jacobi.
# Ejercicio Chapra Ej:6.11
# Sistemas de ecuaciones no lineales
# con método de Newton Raphson para xy
import numpy as np
import sympy as sym
def matrizJacobiano(variables, funciones):
n = len(funciones)
m = len(variables)
# matriz Jacobiano inicia con ceros
Jcb = sym.zeros(n,m)
for i in range(0,n,1):
unafi = sym.sympify(funciones[i])
for j in range(0,m,1):
unavariable = variables[j]
Jcb[i,j] = sym.diff(unafi, unavariable)
return Jcb
# PROGRAMA ----------
# INGRESO
x = sym.Symbol('x')
y = sym.Symbol('y')
f1 = x**2 + x*y - 10
f2 = y + 3*x*(y**2)-57
x0 = 1.5
y0 = 3.5
tolera = 0.0001
# PROCEDIMIENTO
funciones = [f1,f2]
variables = [x,y]
n = len(funciones)
m = len(variables)
Jxy = matrizJacobiano(variables, funciones)
# valores iniciales
xi = x0
yi = y0
# tramo inicial, mayor que tolerancia
itera = 0
tramo = tolera*2
while (tramo>tolera):
J = Jxy.subs([(x,xi),(y,yi)])
# determinante de J
Jn = np.array(J,dtype=float)
determinante = np.linalg.det(Jn)
# iteraciones
f1i = f1.subs([(x,xi),(y,yi)])
f2i = f2.subs([(x,xi),(y,yi)])
numerador1 = f1i*Jn[n-1,m-1]-f2i*Jn[0,m-1]
xi1 = xi - numerador1/determinante
numerador2 = f2i*Jn[0,0]-f1i*Jn[n-1,0]
yi1 = yi -numerador2/determinante
tramo = np.max(np.abs([xi1-xi,yi1-yi]))
xi = xi1
yi = yi1
itera = itera +1
print('iteración: ',itera)
print('Jacobiano con puntos iniciales: ')
sym.pprint(J)
print('determinante: ', determinante)
print('puntos xi,yi:',xi,yi)
print('error:',tramo)
# SALIDA
print('Resultado: ')
print(xi,yi)