{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Systems of Nonlinear Equations\n",
"## CH EN 2450 - Numerical Methods\n",
"**Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah**\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 1\n",
"\n",
"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\n",
"\\begin{equation}\n",
"y = 4x - 0.5 x^3\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"y = \\sin(x)e^{-x}\n",
"\\end{equation}\n",
"\n",
"The first step is to write these in residual form\n",
"\\begin{equation}\n",
"f_1 = y - 4x + 0.5 x^3,\\\\\n",
"f_2 = y - \\sin(x)e^{-x}\n",
"\\end{equation}\n"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy import cos, sin, pi, exp\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import matplotlib.pyplot as plt\n",
"from scipy.optimize import fsolve"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"y1 = lambda x: 4 * x - 0.5 * x**3\n",
"y2 = lambda x: sin(x)*exp(-x)\n",
"x = np.linspace(-3.5,4,100)\n",
"plt.ylim(-8,6)\n",
"plt.plot(x,y1(x), 'k')\n",
"plt.plot(x,y2(x), 'r')\n",
"plt.grid()\n",
"plt.savefig('example1.pdf')"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"def newton_solver(F, J, x, tol):\n",
" F_value = F(x)\n",
" err = np.linalg.norm(F_value, ord=2) # l2 norm of vector\n",
" err = tol + 100\n",
" niter = 0\n",
" while abs(err) > tol and niter < 100:\n",
" J_value = J(x)\n",
" delta = np.linalg.solve(J_value, -F_value)\n",
" x = x + delta\n",
" F_value = F(x)\n",
" err = np.linalg.norm(F_value, ord=2)\n",
" niter += 1\n",
"\n",
" # Here, either a solution is found, or too many iterations\n",
" if abs(err) > tol:\n",
" niter = -1\n",
" return x, niter, err"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"def F(xval):\n",
" x = xval[0]\n",
" y = xval[1]\n",
" f1 = 0.5 * x**3 + y - 4*x\n",
" f2 = y - sin(x)*exp(-x)\n",
" return np.array([f1,f2])\n",
"\n",
"\n",
"def J(xval):\n",
" x = xval[0]\n",
" y = xval[1]\n",
" return np.array([[1.5 * x**2 - 4, 1],\n",
" [-cos(x)*exp(-x)+sin(x)*exp(-x), 1]])\n",
"\n",
"def Jdiff(F,xval):\n",
" n = len(xval)\n",
" J = np.zeros((n,n))\n",
" col = 0\n",
" for x0 in xval:\n",
" delta = np.zeros(n)\n",
" delta[col] = 1e-4 * x0 + 1e-12\n",
" diff = (F(xval + delta) - F(xval))/delta[col]\n",
" J[:,col] = diff\n",
" col += 1\n",
" return J"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4 [-1.46110592 -4.28481694]\n",
"Error Norm = 1.523439979143842e-11\n"
]
}
],
"source": [
"tol = 1e-8\n",
"xguess = np.array([-2,-4])\n",
"x, n, err = newton_solver(F, J, xguess, tol)\n",
"print (n, x)\n",
"print ('Error Norm =',err)"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-1.46110592, -4.28481694])"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fsolve(F,(-2,-4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 2\n",
"Find the roots of the following system of equations\n",
"\\begin{equation}\n",
"x^2 + y^2 = 1, \\\\\n",
"y = x^3 - x + 1\n",
"\\end{equation}\n",
"First we assign $x_1 \\equiv x$ and $x_2 \\equiv y$ and rewrite the system in residual form\n",
"\\begin{equation}\n",
"f_1(x_1,x_2) = x_1^2 + x_2^2 - 1, \\\\\n",
"f_2(x_1,x_2) = x_1^3 - x_1 - x_2 + 1\n",
"\\end{equation}\n"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-1,1)\n",
"y1 = lambda x: x**3 - x + 1\n",
"y2 = lambda x: np.sqrt(1 - x**2)\n",
"plt.plot(x,y1(x), 'k')\n",
"plt.plot(x,y2(x), 'r')\n",
"plt.grid()"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"def F(xval):\n",
" x1 = xval[0]\n",
" x2 = xval[1]\n",
" f1 = x1**2 + x2**2 - 1\n",
" f2 = x1**3 - x1 + 1 - x2\n",
" return np.array([f1,f2])\n",
"\n",
"\n",
"def J(xval):\n",
" x1 = xval[0]\n",
" x2 = xval[1]\n",
" return np.array([[2.0*x1, 2.0*x2],\n",
" [3.0*x1**2 -1, -1]])"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6 [0.74419654 0.66796071]\n",
"Error Norm = 4.965068306494546e-16\n"
]
}
],
"source": [
"tol = 1e-8\n",
"xguess = np.array([0.5,0.5])\n",
"x, n, err = newton_solver(F, J, xguess, tol)\n",
"print (n, x)\n",
"print ('Error Norm =',err)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.74419654, 0.66796071])"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fsolve(F,(0.5,0.5))"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"CSS style adapted from https://github.com/barbagroup/CFDPython. Copyright (c) Barba group\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import urllib\n",
"import requests\n",
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = requests.get(\"https://raw.githubusercontent.com/saadtony/NumericalMethods/master/styles/custom.css\")\n",
" return HTML(styles.text)\n",
"css_styling()\n"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [default]",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
},
"toc": {
"colors": {
"hover_highlight": "#DAA520",
"navigate_num": "#000000",
"navigate_text": "#333333",
"running_highlight": "#FF0000",
"selected_highlight": "#FFD700",
"sidebar_border": "#EEEEEE",
"wrapper_background": "#FFFFFF"
},
"moveMenuLeft": true,
"nav_menu": {
"height": "12px",
"width": "252px"
},
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 4,
"toc_cell": false,
"toc_section_display": "block",
"toc_window_display": false,
"widenNotebook": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}