{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Root finding by homotopy method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We find solution of $f(x) = 0$ by homotopy method applied to\n",
"$$\n",
"h(t,x) = t f(x) + (1-t)(f(x) - f(x_0))\n",
"$$\n",
"We solve the ODE\n",
"$$\n",
"x'(t) = -[h_x(t,x)]^{-1} h_t(t,x), \\qquad x(0) = x_0\n",
"$$\n",
"Let $x = (\\xi_1, \\xi_2) \\in R^2$. Consider\n",
"$$\n",
"f(x) = \\begin{bmatrix}\n",
"\\xi_1^2 - 3 \\xi_2^2 + 3 \\\\\n",
"\\xi_1 \\xi_2 + 6 \\end{bmatrix}\n",
"$$\n",
"Let $x_0 = (1,1)$. Then\n",
"$$\n",
"h_x = f'(x) = \\begin{bmatrix}\n",
"2\\xi_1 & - 6 \\xi_2 \\\\\n",
"\\xi_2 & \\xi_1 \\end{bmatrix}, \\qquad h_t = f(x_0) = \\begin{bmatrix}\n",
"1 \\\\\n",
"7 \\end{bmatrix}\n",
"$$\n",
"$$\n",
"h_x^{-1} = \\frac{1}{\\Delta} \\begin{bmatrix}\n",
"\\xi_1 & 6 \\xi_2 \\\\\n",
"-\\xi_2 & 2\\xi_1 \\end{bmatrix}, \\qquad \\Delta = 2\\xi_1^2 + 6 \\xi_2^2\n",
"$$\n",
"The ODE is given by\n",
"$$\n",
"\\frac{d}{dt}\\begin{bmatrix}\n",
"\\xi_1 \\\\ \\xi_2 \\end{bmatrix} = -h_x^{-1} h_t =\n",
"-\\frac{1}{\\Delta} \\begin{bmatrix}\n",
"\\xi_1 + 42 \\xi_2 \\\\\n",
"-\\xi_2 + 14 \\xi_1 \\end{bmatrix}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import numpy as np\n",
"from scipy.integrate import odeint\n",
"from scipy.linalg import solve, norm\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the function whose roots are required."
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" y = np.zeros(2)\n",
" y[0] = x[0]**2 - 3*x[1]**2 + 3\n",
" y[1] = x[0]*x[1] + 6\n",
" return y\n",
"\n",
"def df(x):\n",
" y = np.array([[2*x[0],-6*x[1]],\n",
" [x[1], x[0]]])\n",
" return y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the rhs of the ODE obtained from homotopy method."
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [],
"source": [
"def F(x,t):\n",
" y = np.zeros(2)\n",
" delta = 2*x[0]**2 + 6*x[1]**2\n",
" y[0] = -(x[0] + 42*x[1])/delta\n",
" y[1] = -(-x[1] + 14*x[0])/delta\n",
" return y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve the ode with a relaxed error tolerance."
]
},
{
"cell_type": "code",
"execution_count": 96,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x = [-2.99046055 2.00603356]\n",
"f(x) = [-0.12965756 0.00103578]\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x0 = np.array([1.0,1.0])\n",
"t = np.linspace(0,1,100)\n",
"x = odeint(F,x0,t,rtol=0.1)\n",
"\n",
"# plot results\n",
"plt.plot(t,x)\n",
"plt.xlabel('t')\n",
"plt.ylabel('x(t)')\n",
"plt.grid(True)\n",
"\n",
"# Final solution\n",
"xf = np.array([x[-1,0],x[-1,1]])\n",
"print('x =',xf)\n",
"print('f(x) =',f(xf))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can improve the final solution by applying Newton-Raphson method."
]
},
{
"cell_type": "code",
"execution_count": 97,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" -2.9904605534e+00 2.0060335560e+00 1.2966169984e-01\n",
" -2.9999822220e+00 1.9999926789e+00 6.0518131105e-05\n",
" -3.0000000000e+00 2.0000000000e+00 2.0260107090e-10\n",
" -3.0000000000e+00 2.0000000000e+00 0.0000000000e+00\n",
" -3.0000000000e+00 2.0000000000e+00 0.0000000000e+00\n"
]
}
],
"source": [
"N = 10\n",
"eps = 1.0e-13\n",
"\n",
"print('%18.10e %18.10e %18.10e' % (xf[0], xf[1], norm(f(xf))))\n",
"for i in range(N):\n",
" J = df(xf)\n",
" dx = solve(J,-f(xf))\n",
" xf = xf + dx\n",
" print('%18.10e %18.10e %18.10e' % (xf[0], xf[1], norm(f(xf))))\n",
" if norm(dx) < eps*norm(xf):\n",
" break"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}