How to solve non-linear equations using python

Question

I have the following code:

#!/usr/bin/env python
from scipy.optimize import fsolve
import math

h = 6.634e-27
k = 1.38e-16

freq1 = 88633.9360e6
freq2 = 88631.8473e6
freq3 = 88630.4157e6

def J(freq,T):
     return (h*freq/k)/(math.exp(h*freq/(k*T))-1)

def equations(x,y,z,w,a,b,c,d):
     f1 = a*(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
     f2 = b*(J(freq3,w)-J(freq3,2.73))*(1-math.exp(-b*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
     f3 = c*(J(freq3,w)-J(freq3,2.73))*(1-math.exp(-b*z))-(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))
     f4 = d*(J((freq3+freq1)/2,(y+w)/2)-J((freq3+freq1)/2,2.73))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
     return (f1,f2,f3,f4)

So, I have defined the equations in the above code. The system of 4-non-linear equations includes 4 variables, a->d, that are predetermined and 4 unknowns, x,y,z and w. I wish to somehow define a->d and feed them into fsolve, thereby creating unique solutions for x,y,z and w. Is this possible?

Answer 1

I'm not getting the x is not defined error you get, but that might be something going wrong on the command line. This can easily be done in a script however with no changes.

Scipy's fsolve takes a callable function as input. So if you want to solve f1 for example, you'd do something like this:

def f1(x, y, z, a):
    return a*(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))

root = fsolve(f1, x0=1, args=(1, 1, 1))

This solves f1(x, y, z, a) = 0 for x with y, z and a as additional arguments (being 1 in this case). So depending on which variable you want to solve for, that should appear first in the argument order e.g. solving for a would require f1(a, x, y, z). The x0 is the starting estimate for the root.

This then returns an ndarray with the root.

For more information check out http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html

Update: Looking at the answer to How to solve a pair of nonlinear equations using Python?

To solve for all of x, y, z, and w you need to pass them as a single argument. So you end up with something like this:

h = 6.634e-27
k = 1.38e-16

freq1 = 88633.9360e6
freq2 = 88631.8473e6
freq3 = 88630.4157e6

def J(freq,T):
    return (h*freq/k)/(math.exp(h*freq/(k*T))-1)

def equations(p, a, b, c, d):
    x, y, z, w = p
    f1 = a*(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
    f2 = b*(J(freq3,w)-J(freq3,2.73))*(1-math.exp(-b*z))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
    f3 = c*(J(freq3,w)-J(freq3,2.73))*(1-math.exp(-b*z))-(J(freq1,y)-J(freq1,2.73))*(1-math.exp(-a*z))
    f4 = d*(J((freq3+freq1)/2,(y+w)/2)-J((freq3+freq1)/2,2.73))-(J(freq2,x)-J(freq2,2.73))*(1-math.exp(-z))
    return (f1,f2,f3,f4)

x, y, z, w = fsolve(equations, x0=(1, 1, 1, 1), args=(1, 1, 1, 1), xtol=1e-13) # x0=(x, y, z, w) args=(a, b, c, d)
print x, y, z, w # prints [2.73, 2.73, <some value>, 2.73]
print equations((x, y, z, w), 1, 1, 1, 1) # prints a tuple with values of order 1e-13 to 1e-16 (so we can assume they are indeed the roots)