In [1]:
# Solving nonlinear systems of equations

# To solve a system of nonlinear equations, we will use 'fsolve()' which requires
# the `scipy.optimize' module.
import numpy as np
from scipy.optimize import fsolve
In [2]:
# First, we will solve a single nonlinear equations.  We start by defining
# our equation as a function.
def f(x):
    return 3*x**3 - 2*x**2 + x - 7
In [4]:
# Once the function is defined, we can evaluate it for different values of x
# by:
f(1)
Out[4]:
-5
In [7]:
# We can now use fsolve().  The first arguement is the function (i.e. the equation
# that we want to solve, assumed to be equal to zero) and the second argument
# is an initial guess at the solution.
x = fsolve(f, 1)
In [8]:
# The output of fsolve() is an array.
x
Out[8]:
array([1.49175209])
In [9]:
# We can index the array to get the final solution.
x[0]
Out[9]:
1.4917520892393685
In [10]:
# Let's check that the solution works.
f(x[0])
Out[10]:
-1.2434497875801753e-14
In [11]:
# Here's another example: Find y such that tan(e^(-2y)) = 1/y.
def f(y):
    return np.tan(np.exp(-2*y)) - 1/y
In [14]:
# Here's the solution to f(y) = 0.
y = fsolve(f, -0.1)
y[0]
Out[14]:
-0.31441845309470834
In [15]:
# Test the solution...
f(y[0])
Out[15]:
-2.220446049250313e-15
In [17]:
# Note that there is another way to use fsolve such that it outputs just the 
# numerical solution (rather than inserting it into an array).
y, = fsolve(f, -0.1)
print('y =', y)
print('f(y) =', f(y))

y = -0.31441845309470834
f(y) = -2.220446049250313e-15
In [18]:
# We can also solve a system of nonlinear equations.  Below we define a system
# of three equations with three unknown variables.
def equations(vars):
    x, y, z = vars
    eq1 = 3*x*y + y - z -10
    eq2 = x + x**2*y + z -12
    eq3 = x -y -z + 2
    return [eq1, eq2, eq3]
In [19]:
# To evaluate the three equations at particular values of x, y, and z, we can
# enter:
equations([1, 0, -3])
Out[19]:
[-7, -14, 6]
In [20]:
# Here is the fsolve() statement with an initial guess at the solution
x, y, z =  fsolve(equations, (1, 1, 1))
In [24]:
# Using this 'comma' notation, the get the three solutions directly.
print('x =', x, 'y =', y, 'z =', z)

x = 2.100953871762653 y = 1.6983245687167783 z = 2.4026293030458747
In [32]:
# Here, we check that all three of the equations actually evaluate to something
# close to zero.  Note that, the 'SymPy' module is needed to use the 'N()' function. 
from sympy import *
print('eq1 =', N(equations([x, y, z])[0], 3), 'eq2 =', N(equations([x, y, z])[1], 3),\
      'eq3 =', N(equations([x, y, z])[2], 3))

eq1 = 1.36e-10 eq2 = 2.76e-10 eq3 = 0