0

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?

1 Answer 1

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)
Sign up to request clarification or add additional context in comments.

2 Comments

This may work, but doesn't solve the problem at all. I require the solution to x,y,z, and w. The values of a,b,c,d are predetermined.
I updated it for multiple variables. The solution for z varies wildly depending on the initial guess x0. I hope this is what you are looking for, otherwise I'm sorry for causing confusion.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.