Radical Edward Radical Edward - 2 months ago 12
Python Question

Solving sets of non-linear equations as an array

I'm trying to solve for the intersection of the two equations:

y=Rx^1.75
and
y=ax^2+bx+c
for all rows in my dataframe (about 100K rows). Each value of
R,a,b,c
is different for each row. I can solve them one by one by iterating through the dataframe and calling
fsolve()
for each row (as done below), but I'm wondering if there is a better way to do this.

My question is: is it possible to turn this into an array calculation, that is, to solve all the rows at once? Any ideas on how to get this calculation done faster would be really helpful.

Here is an example dataframe with coefficients

R a b c
0 0.5 -0.01 -0.50 32.42
1 0.6 0.00 0.07 14.12
2 0.7 -0.01 -0.50 32.42


And here is the working example code that I'm using to test methods:

import numpy as np
import pandas as pd
from scipy.optimize import *

# The fSolve function
def myFunction(zGuess,*Params):
# Get the coefficients
R,a,b,c = Params
# Get the initial guess
x,y = zGuess

F = np.empty((2))
F[0] = R*x**1.75-y
F[1] = a*x**2+b*x+c-y
return F

# Example Dataframe that is 10K rows of different coefficients
df = pd.DataFrame({"R":[0.500, 0.600,0.700],
"a":[-0.01, 0.000,-0.01],
"b":[-0.50, 0.070,-0.50],
"c":[32.42, 14.12,32.42]})
# Initial guess
zGuess = np.array([50,50])

# Make a place to store the answers
df["x"] = None
df["y"] = None

# Loop through the rows?
for index, coeffs in df.iterrows():
# Get the coefficients
Params = (coeffs["R"],coeffs["a"],coeffs["b"],coeffs["c"])
# fSolve
z = fsolve(myFunction,zGuess,args=Params)
# Set the answers
df.loc[index,"x"] = z[0]
df.loc[index,"y"] = z[1]

print df


============================================

Solution (who's answer is faster):



I got two answers below that both gave mathematically correct answers. So at this point, it's all who's calculation is faster! The test dataframe will be 3K rows.

Answer #1 (Newtons Method)

# Solution 1
import numpy as np
import pandas as pd
Count = 1000
df = pd.DataFrame({"R":[0.500, 0.600,0.700]*Count,
"a":[-0.01, 0.000,-0.01]*Count,
"b":[-0.50, 0.070,-0.50]*Count,
"c":[32.42, 14.12,32.42]*Count})

from datetime import datetime
t_start = datetime.now()
#---------------------------------
InitialGuess = 50.0
Iterations = 20
x = np.full(df["a"].shape, InitialGuess)

for i in range(Iterations):
x = x - (-df["R"]*x**1.75 + df["a"]*x**2 + df["b"]*x + df["c"])/(-1.75*df["R"]*x**0.75 + 2*df["a"]*x + df["b"])

df["x"] = x
df["y"] = df["R"]*x**1.75
df["x Error"] = df["a"]*x**2 + df["b"]*x + df["c"] - df["R"]*x**1.75
#---------------------------------
t_end = datetime.now()
print ('\n\n\nTime spent running this was:')
print(t_end - t_start)
print df


And the time spent was:

Time spent running this was:
0:00:00.015000


Answer #2 (fSolve)

# Solution 2
import numpy as np
import pandas as pd
from scipy.optimize import *
Count = 1000
df = pd.DataFrame({"R":[0.500, 0.600,0.700]*Count,
"a":[-0.01, 0.000,-0.01]*Count,
"b":[-0.50, 0.070,-0.50]*Count,
"c":[32.42, 14.12,32.42]*Count})

from datetime import datetime
t_start = datetime.now()
#---------------------------------
coefs = df.values[:, 0:4]

def mfun(x, *args):
args = np.array(args[0], dtype=np.float64)
return args[:,1] * x**2 + args[:,2] * x + args[:,3] - args[:,0] * x**1.75


nrows = coefs.shape[0]
df["x"] = fsolve(mfun, np.ones(nrows) * 50, args=coefs)
df["y"] = coefs[:, 0] * df["x"]**1.75
#---------------------------------
t_end = datetime.now()
print ('\n\n\nTime spent running this was:')
print(t_end - t_start)
print df


And the time spent was:

Time spent running this was:
0:00:35.786000


Final thoughts:



For this particular case, Newtons method was much faster (I can run 300K rows in
0:00:01.139000
!). Thank you both!

Answer

Maybe you can use Newton's method:

import numpy as np

data = np.array(
    [[0.5, -0.01, -0.50,  32.42],
     [0.6,  0.00,  0.07,  14.12],
     [0.7, -0.01, -0.50,  32.42]])


R, a, b, c = data.T

x = np.full(a.shape, 10.0)
m = 1.0
for i in range(20):
    x = x - m * (-R*x**1.75 + a*x**2 + b*x + c)/(-1.75*R*x**0.75 + 2*a*x + b)

print(a*x**2 + b*x + c - R * x**1.75)

output:

[  0.00000000e+00   1.77635684e-15   3.55271368e-15]

be careful to choose the iteration count and initial value of x.

Comments