I've been working with this equations of Potential Flow
pSI: 0 = 10 * y + 23.8732414638 * (pi + atan(abs((x - 12) / x))) + 23.8732414638 * (pi - atan(abs((y + 12) / x)))-234.882642242
V: 0 = ((10 + 23.8732414638 * x / (x*2 + (y - 12)*2) + 23.8732414638 * x / (x*2 + (y + 12)*2))**2 + (23.8732414638 * (y - 12) / (x*2 + (y - 12)*2) + 23.8732414638 * (y + 12) / (x*2 + (y + 12)*2))*2)*0.5-8
using the next code to solve them:
# Prototype of N-R for a system of two non-linear equations
# evaluating functions of two variables
from math import *
eq1 = raw_input('Enter the equation 1: ')
eq2 = raw_input('Enter the equation 2: ')
x0 = float(input('Enter x: '))
y0 = float(input('Enter y: '))
def f(x,y):
return eval(eq1)
def g(x,y):
return eval(eq2)
x = x0
y = y0
for n in range(1, 8):
a = (f(x + 1e-06, y) - f(x,y)) / 1e-06
b = (f(x, y + 1e-06) - f(x,y)) / 1e-06
c = - f(x,y)
d = (g(x + 1e-06, y) - g(x,y)) / 1e-06
eE = (g(x, y + 1e-06) - g(x,y)) / 1e-06
fF = - g(x,y)
print "f(x, y)= ", eq1
print "g(x, y)= ", eq2
print """x y """
print x, y
print """a b c d e f """
print a, b, c, d, e, fF
print """
a * x + b * y = c
d * x + e * y = f
"""
print a," * x + ",b," * y = ",c
print d," * x + ",eE," * y = ",fF
_Sy = (c - a * fF / d) / (b - a * eE / d)
_Sx = (fF / d) - (eE / d) * _Sy
Ea_X = (_Sx ** 2 + _Sy ** 2)**0.5
x = x + _Sx
y = y + _Sy
print "Sx = ", _Sx
print "Sy = ", _Sy
print "x = ", x
print "y = ", y
print "|X_1 - X_0| = ", Ea_X
And gives me the Zero division Error
Traceback (most recent call last):
File "C:\Documents and Settings\Principal\Mis documentos\Finished_Non_lin_N_R.py", line 38, in <module>
d = (g(x + 1e-06, y) - g(x,y)) / 1e-06
File "C:\Documents and Settings\Principal\Mis documentos\Finished_Non_lin_N_R.py", line 27, in g
return eval(eq2)
File "<string>", line 1, in <module>
ZeroDivisionError: float division by zero
Then I tried with this one that I've found in one solved question:
from scipy.optimize import fsolve
from math import *
def equations(p):
x, y = p
return ( ((10 + 23.8732414638 * x / (x**2 + (y - 12)**2) + 23.8732414638 * x / (x**2 + (y + 12)**2))**2 + (23.8732414638 * (y -12) / (x**2 + (y - 12)**2) + 23.8732414638 * (y + 12) / (x**2 + (y + 12)**2))**2)**0.5-8, 10 * y + 23.8732414638 * (pi + atan(abs((x - 12) / x))) + 23.8732414638 * (pi - atan(abs((y + 12) / x)))-234.882642242)
x, y = fsolve(equations, (0, 18))
print equations((x, y))
And with the next problem:
<module2>:19: RuntimeWarning: divide by zero encountered in long_scalars
<module2>:19: RuntimeWarning: divide by zero encountered in double_scalars
(-1.3374190643844486e-11, -2.8308022592682391e-11)
And the last part of this story is that I have some code that actually works but is too clumsy and a little bit inaccurate, and I would like some suggestions in order to make it simple:
from math import *
U = 10
b = 15
h = 12
cF = 0.5 * U * b / pi
pSI0 = 234.882642242
VP = 12.0
X0 = 0
Y0 = 40
sX = 0.01
sY = 0.01
print cF
def FirstFor(A,B,C,D):
for n in range(1,808):
for m in range(1,4002):
X = A - B*n
Y = C - D*m
pSI = U * Y + cF * (pi + atan(abs((Y - h) / X))) + cF * (pi - atan(abs((Y + h) / X)))
Vaprox = ((10 + 23.8732414638 * X / (X**2 + (Y - 12)**2) + 23.8732414638 * X / (X**2 + (Y + 12)**2))**2 +(23.8732414638 * (Y - 12) / (X**2 + (Y - 12)**2) + 23.8732414638 * (Y + 12) / (X**2 + (Y + 12)**2))**2)**0.5
EA1 = abs(pSI0 - pSI)
EA2 = abs(VP - Vaprox)
if EA1 < 0.05 and EA2 < 0.05:
X1f = X
Y1f = Y
H3 = [X1f,Y1f]
print n,m,X,Y,"GG son",EA1,EA2
print H3
for n in range(1,808):
for m in range(1,4002):
X = H3[0] - B*n*0.001
Y = H3[1] - D*m*0.001
pSI = U * Y + cF * (pi + atan(abs((Y - h) / X))) + cF * (pi - atan(abs((Y + h) / X)))
Vaprox = ((10 + 23.8732414638 * X / (X**2 + (Y - 12)**2) + 23.8732414638 * X / (X**2 + (Y + 12)**2))**2 +(23.8732414638 * (Y - 12) / (X**2 + (Y - 12)**2) + 23.8732414638 * (Y + 12) / (X**2 + (Y + 12)**2))**2)**0.5
EA1 = abs(pSI0 - pSI)
EA2 = abs(VP - Vaprox)
if EA1 < 0.0001 and EA2 < 0.0001:
X2f = X
Y2f = Y
I3 = [X2f,Y2f]
print n,m,X,Y,"GG son",EA1,EA2
print I3
for n in range(1,808):
for m in range(1,4002):
X = H3[0] + B*n*0.001
Y = H3[1] + D*m*0.001
pSI = U * Y + cF * (pi + atan(abs((Y - h) / X))) + cF * (pi - atan(abs((Y + h) / X)))
Vaprox = ((10 + 23.8732414638 * X / (X**2 + (Y - 12)**2) + 23.8732414638 * X / (X**2 + (Y + 12)**2))**2 +(23.8732414638 * (Y - 12) / (X**2 + (Y - 12)**2) + 23.8732414638 * (Y + 12) / (X**2 + (Y + 12)**2))**2)**0.5
EA1 = abs(pSI0 - pSI)
EA2 = abs(VP - Vaprox)
if EA1 < 0.0001 and EA2 < 0.0001:
X2f = X
Y2f = Y
I3 = [X2f,Y2f]
print n,m,X,Y,"GG son",EA1,EA2
print I3
# starting with the FirstFor
FirstFor(X0,sX,Y0,sY)
print "\n This should be good"
raw_input()