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.
Instrucciones 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)
