#!/usr/bin/env python
# coding: utf-8
# # Systems of Nonlinear Equations
# ## CH EN 2450 - Numerical Methods
# **Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah**
#
# # Example 1
#
# A system of nonlinear equations consists of several nonlinear functions - as many as there are unknowns. Solving a system of nonlinear equations means funding those points where the functions intersect each other. Consider for example the following system of equations
# \begin{equation}
# y = 4x - 0.5 x^3
# \end{equation}
# \begin{equation}
# y = \sin(x)e^{-x}
# \end{equation}
#
# The first step is to write these in residual form
# \begin{equation}
# f_1 = y - 4x + 0.5 x^3,\\
# f_2 = y - \sin(x)e^{-x}
# \end{equation}
#
# In[21]:
import numpy as np
from numpy import cos, sin, pi, exp
get_ipython().run_line_magic('matplotlib', 'inline')
get_ipython().run_line_magic('config', "InlineBackend.figure_format = 'svg'")
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
# In[22]:
y1 = lambda x: 4 * x - 0.5 * x**3
y2 = lambda x: sin(x)*exp(-x)
x = np.linspace(-3.5,4,100)
plt.ylim(-8,6)
plt.plot(x,y1(x), 'k')
plt.plot(x,y2(x), 'r')
plt.grid()
plt.savefig('example1.pdf')
# In[26]:
def newton_solver(F, J, x, tol):
F_value = F(x)
err = np.linalg.norm(F_value, ord=2) # l2 norm of vector
err = tol + 100
niter = 0
while abs(err) > tol and niter < 100:
J_value = J(x)
delta = np.linalg.solve(J_value, -F_value)
x = x + delta
F_value = F(x)
err = np.linalg.norm(F_value, ord=2)
niter += 1
# Here, either a solution is found, or too many iterations
if abs(err) > tol:
niter = -1
return x, niter, err
# In[27]:
def F(xval):
x = xval[0]
y = xval[1]
f1 = 0.5 * x**3 + y - 4*x
f2 = y - sin(x)*exp(-x)
return np.array([f1,f2])
def J(xval):
x = xval[0]
y = xval[1]
return np.array([[1.5 * x**2 - 4, 1],
[-cos(x)*exp(-x)+sin(x)*exp(-x), 1]])
def Jdiff(F,xval):
n = len(xval)
J = np.zeros((n,n))
col = 0
for x0 in xval:
delta = np.zeros(n)
delta[col] = 1e-4 * x0 + 1e-12
diff = (F(xval + delta) - F(xval))/delta[col]
J[:,col] = diff
col += 1
return J
# In[28]:
tol = 1e-8
xguess = np.array([-2,-4])
x, n, err = newton_solver(F, J, xguess, tol)
print (n, x)
print ('Error Norm =',err)
# In[29]:
fsolve(F,(-2,-4))
# # Example 2
# Find the roots of the following system of equations
# \begin{equation}
# x^2 + y^2 = 1, \\
# y = x^3 - x + 1
# \end{equation}
# First we assign $x_1 \equiv x$ and $x_2 \equiv y$ and rewrite the system in residual form
# \begin{equation}
# f_1(x_1,x_2) = x_1^2 + x_2^2 - 1, \\
# f_2(x_1,x_2) = x_1^3 - x_1 - x_2 + 1
# \end{equation}
#
# In[39]:
x = np.linspace(-1,1)
y1 = lambda x: x**3 - x + 1
y2 = lambda x: np.sqrt(1 - x**2)
plt.plot(x,y1(x), 'k')
plt.plot(x,y2(x), 'r')
plt.grid()
# In[33]:
def F(xval):
x1 = xval[0]
x2 = xval[1]
f1 = x1**2 + x2**2 - 1
f2 = x1**3 - x1 + 1 - x2
return np.array([f1,f2])
def J(xval):
x1 = xval[0]
x2 = xval[1]
return np.array([[2.0*x1, 2.0*x2],
[3.0*x1**2 -1, -1]])
# In[31]:
tol = 1e-8
xguess = np.array([0.5,0.5])
x, n, err = newton_solver(F, J, xguess, tol)
print (n, x)
print ('Error Norm =',err)
# In[32]:
fsolve(F,(0.5,0.5))
# In[40]:
import urllib
import requests
from IPython.core.display import HTML
def css_styling():
styles = requests.get("https://raw.githubusercontent.com/saadtony/NumericalMethods/master/styles/custom.css")
return HTML(styles.text)
css_styling()