{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Lagrange Multiplier",
   "id": "9f8d2321e74ad9dc"
  },
  {
   "cell_type": "code",
   "id": "initial_id",
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2026-04-28T11:35:33.360292Z",
     "start_time": "2026-04-28T11:35:32.570461Z"
    }
   },
   "source": [
    "import numpy as np\n",
    "import sympy as sp"
   ],
   "outputs": [],
   "execution_count": 1
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Example",
   "id": "7c20770541717344"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-04-28T11:50:46.463985Z",
     "start_time": "2026-04-28T11:50:46.449587Z"
    }
   },
   "cell_type": "code",
   "source": [
    "x, y, la= sp.symbols(\"x y la\", real = True)\n",
    "\n",
    "f = (x - 3)**2 + (y - 2)**2\n",
    "f_num = sp.lambdify((x, y), f)\n",
    "phi = (x / 3)**2 + (y / 2)**2 - 1\n",
    "phi_num = sp.lambdify((x, y), phi)\n",
    "L = f - la * phi\n",
    "L_num = sp.lambdify((x, y, la), L)\n",
    "\n",
    "L_grad = sp.Matrix([sp.diff(L, x), sp.diff(L, y), sp.diff(L, la)])\n",
    "L_grad_num = sp.lambdify((x, y, la), list(L_grad)) # use list here to get a nicer output as flat array\n",
    "\n",
    "H = sp.hessian(L, (x, y, la))\n",
    "H_num = sp.lambdify((x, y, la), H)"
   ],
   "id": "ab90b1b2981efcbd",
   "outputs": [],
   "execution_count": 35
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### Newton's Method",
   "id": "7ae85714f38e30d2"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-04-28T11:53:30.216830Z",
     "start_time": "2026-04-28T11:53:30.210142Z"
    }
   },
   "cell_type": "code",
   "source": [
    "sol_old = np.array([2, 1, -2]) # for min\n",
    "sol_new = np.array([2, 1, -2]) # for min\n",
    "\n",
    "sol_old = np.array([-2, -1, -2]) # for max\n",
    "sol_new = np.array([-2, -1, -2]) # for max\n",
    "solutions = [sol_old]\n",
    "\n",
    "for i in range(20):\n",
    "    g_old = np.array(L_grad_num(*sol_old))\n",
    "    H_old = np.array(H_num(*sol_old))\n",
    "    sol_new = sol_old - np.linalg.solve(H_old, g_old) # better to avoid inverse and ask numpy to solve a system of linear equations directly\n",
    "    solutions.append(sol_new)\n",
    "    sol_old = sol_new\n",
    "    \n",
    "sol_new"
   ],
   "id": "b7048cfe7a2a52ed",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.88143787, -0.55670291, 18.37032177])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 39
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-04-28T11:54:00.205010Z",
     "start_time": "2026-04-28T11:54:00.199887Z"
    }
   },
   "cell_type": "code",
   "source": "np.sqrt(f_num(sol_new[0], sol_new[1]))",
   "id": "8541e3b5307bbd4d",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.413114781967053"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 40
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-04-28T11:51:51.453305Z",
     "start_time": "2026-04-28T11:51:51.447631Z"
    }
   },
   "cell_type": "code",
   "source": "phi_num(sol_new[0], sol_new[1])",
   "id": "80bc3febde1805c",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 38
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}