求解2个非线性方程组

Question

我一直在研究这个势流方程

pSI:0=10*y+23.8732414638*（pi+atan（绝对值（（x-12）/x）））+23.8732414638*（pi-atan（绝对值（（y+12）/x））-234.88264242

V:0=（（10+23.8732414638*x/（x*2+（y-12）*2）+23.8732414638*x/（x*2+（y+12）（y+12*2）2））**2+（23.8732414638*（y-12）/（x*2+（y-12）*2）+23.8732414638*（y+12）/（x*2+（y+12）*2））））*2）[y+12）*2））*2）2）2）*0.5-8）0.5-8）0在

使用下一个代码来解决这些问题：

# 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

给我零除法误差

^{pr2}$

然后我试着在一个已解决的问题中找到了这个：

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

接下来的问题是：

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

这个故事的最后一部分是，我有一些实际可行的代码 但是太笨拙了，而且有点不准确，为了使其简单化，我想提出一些建议：

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()