Chapter-20.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "from sympy import *\n",
    "init_printing(use_latex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 20 - Vector Fields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20.1.1 - Coordinate Systems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.vector import *\n",
    "N = CoordSys3D(\"N\")\n",
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.i, N.j, N.k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = symbols(\"t\")\n",
    "v1 = (t**2 + 1) * N.i\n",
    "v2 = (t**2 + 1) * N.i - 2 * N.j + N.k\n",
    "display(v1, v2)\n",
    "print(type(v1), type(v2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i, j, k = N.base_vectors()\n",
    "v3 = i + 2 * j + 3 * k\n",
    "v3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.x, N.y, N.z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, z = N.base_scalars()\n",
    "f = x**2 * y + z\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "O = CoordSys3D(\"\", vector_names=[\"e_1\", \"e_2\", \"e_3\"], variable_names=[\"a\", \"b\", \"c\"])\n",
    "a, b, c = O.base_scalars()\n",
    "e1, e2, e3 = O.base_vectors()\n",
    "v = (a + b) * e1 + (b + c) * e2 + (a + c) * e3\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = CoordSys3D(\"M\", transformation=\"cylindrical\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = CoordSys3D('A') # initial Cartesian CS\n",
    "B = A.locate_new(\"B\", position=A.i + A.j + A.k) # intermediate Cartesian CS\n",
    "C = B.create_new(\"C\", transformation=\"cylindrical\")\n",
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = CoordSys3D(\"N\")\n",
    "M = N.create_new(\"M\", transformation=\"cylindrical\")\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M.base_scalars()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(M.transformation_to_parent(), M.transformation_from_parent())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i, j, k = N.base_vectors()\n",
    "# write v in N\n",
    "v = 2 * i + 3 * j + 4 * k\n",
    "# express v in M\n",
    "express(v, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# substitution dictionary\n",
    "d = {key: val for key, val in zip(N.base_scalars(), v.components.values())}\n",
    "# values of the componenents in the new reference system\n",
    "comps = list((c.subs(d) for c in M.transformation_from_parent()))\n",
    "# combining base vectors with the new components\n",
    "new_comps = list((b * c for b, c in zip(M.base_vectors(), comps)))\n",
    "# final result\n",
    "new_v = VectorAdd(*new_comps)\n",
    "new_v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = v.to_matrix(N)\n",
    "display(A)\n",
    "print(type(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_to_vector(A, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20.1.2 - Vector Operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = CoordSys3D(\"N\")\n",
    "x, y, z = N.base_scalars()\n",
    "i, j, k = N.base_vectors()\n",
    "a, b, c = symbols(\"a:c\")\n",
    "v1 = x * y * i + y * z * j + z * x * k\n",
    "v2 = a * i + b * j + c * k\n",
    "display(v1, v2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1 + v2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2 * v2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(v1.dot(v2), v1 ^ v2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "field = x**2 * y + z\n",
    "gradient(field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "divergence(v1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "curl(v1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v3 = v1 * x\n",
    "laplacian(v3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "directional_derivative(field, v1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = Del()\n",
    "(d & v1).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(v1 & d)(field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = x * i + (x + y) * j + k\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.integrate(x).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.diff(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 20.2 - Vector Integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.abc import t, theta, phi, x, y, z\n",
    "from sympy.geometry import Point\n",
    "from sympy.vector import vector_integrate, ParametricRegion, ImplicitRegion\n",
    "r, H, R = symbols(\"r, H, R\", real=True, positive=True)\n",
    "circle1 = ParametricRegion(\n",
    "    (r * cos(theta) - 2, r * sin(theta) + 3), (theta, 0, 2 * pi))\n",
    "circle2 = ImplicitRegion((x, y), (x + 2)**2 + (y - 3)**2 - r**2)\n",
    "circle3 = Circle(Point(-2, 3), r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "disc = ParametricRegion(\n",
    "    (r * cos(theta), r * sin(theta)), (r, 0, 2), (theta, 0, 2 * pi))\n",
    "annulus = ParametricRegion(\n",
    "    (r * cos(theta), r * sin(theta)), (r, 1, 2), (theta, 0, 2 * pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cylinder_surface = ParametricRegion((r * cos(theta), r * sin(theta), z),\n",
    "    (z, 0, H), (theta, 0, 2 * pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cylinder_vol = ParametricRegion((r * cos(theta), r * sin(theta), z),\n",
    "    (z, 0, H), (r, 0, R), (theta, 0, 2 * pi))\n",
    "annular_cross_section_vol = ParametricRegion((r * cos(theta), r * sin(theta), z),\n",
    "    (z, 0, H), (r, 0, R), (theta, 0, 2 * pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "segment = Segment(Point(4, -1, 9), Point(1, 5, 7))\n",
    "triangle = Polygon((0, 0), (1, 0), (1, 1))\n",
    "box = ParametricRegion((x, y, z), (x, -1, 1), (y, -2, 2), (z, -3, 3))\n",
    "cone = ParametricRegion(\n",
    "    (r * cos(theta), r * sin(theta), r), (r, 0, 3), (theta, 0, 2 * pi))\n",
    "parabola = ImplicitRegion((x, y), (y - 1)**2 - 4*(x + 6))\n",
    "ellipse = ImplicitRegion((x, y), (x**2 / 4 + y**2 / 16 - 1))\n",
    "sphere1 = ImplicitRegion((x, y, z), x**2 + y**2 + z**2 - 2 * x)\n",
    "sphere2 = ParametricRegion(\n",
    "    (r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)),\n",
    "    (phi, 0, pi/2), (theta, 0, 2*pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere1.rational_parametrization()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_curve = vector_integrate(1, circle1)\n",
    "area_disc = vector_integrate(1, disc)\n",
    "volume_cyl = vector_integrate(1, cylinder_vol)\n",
    "display(length_curve, area_disc, volume_cyl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20.2.1 - Line Integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from sympy.vector import *\n",
    "t = symbols(\"t\")\n",
    "C = CoordSys3D(\"C\")\n",
    "x, y, z = C.base_scalars()\n",
    "half_circle = ParametricRegion((4 * cos(t), 4 * sin(t)), (t, 3*pi/2, pi/2))\n",
    "vector_integrate(x*y**4, half_circle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "helix = ParametricRegion((cos(t), sin(t), 3*t), (t, 0, 4*pi))\n",
    "vector_integrate(x*y*z, helix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20.2.2 - Line Integrals of Vector Fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = CoordSys3D(\"C\")\n",
    "x, y, z = C.base_scalars()\n",
    "i, j, k = C.base_vectors()\n",
    "\n",
    "F = x*z*i -y*z*k\n",
    "curve = Segment3D((-1, 2, 0), (3, 0, 1))\n",
    "vector_integrate(F, curve)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20.2.3 - Surface Integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = CoordSys3D(\"C\")\n",
    "i, j, k = C.base_vectors()\n",
    "x, y, z = C.base_scalars()\n",
    "\n",
    "rs, phi, theta = symbols(\"r, phi, theta\")\n",
    "r = sqrt(3)\n",
    "\n",
    "cylinder = ParametricRegion(\n",
    "    (r * cos(theta), r * sin(theta), z),\n",
    "    (z, 0, 4 - r * sin(theta)), (theta, 0, 2*pi))\n",
    "top_plane = ParametricRegion(\n",
    "    (rs * cos(theta), rs * sin(theta), 4 - rs * sin(theta)),\n",
    "    (theta, 0, 2*pi), (rs, 0, r))\n",
    "bottom_plane = ParametricRegion(\n",
    "    (rs * cos(theta), rs * sin(theta), 0),\n",
    "    (theta, 0, 2*pi), (rs, 0, r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from spb import *\n",
    "\n",
    "sliced_cylinder = Piecewise(\n",
    "    (z, z < 4 - r * sin(theta)),\n",
    "    (4 - r * sin(theta), True)\n",
    ")\n",
    "display(sliced_cylinder)\n",
    "\n",
    "options = dict(wireframe=True, wf_n1=20)\n",
    "(tp_fx, tp_fy, tp_fz), tp_r1, tp_r2 = top_plane.args\n",
    "(bp_fx, bp_fy, bp_fz), bp_r1, bp_r2 = bottom_plane.args\n",
    "\n",
    "graphics(\n",
    "    surface_parametric(r * cos(theta), r * sin(theta), sliced_cylinder,\n",
    "        (theta, 0, 2*pi), (z, 0, 8), **options),\n",
    "    surface_parametric(tp_fx, tp_fy, tp_fz, tp_r1, tp_r2, **options),\n",
    "    surface_parametric(bp_fx, bp_fy, bp_fz, bp_r1, bp_r2, **options),\n",
    "    backend=PB, legend=False\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "field = y + z\n",
    "r1 = vector_integrate(field, cylinder)\n",
    "r2 = vector_integrate(field, top_plane)\n",
    "r3 = vector_integrate(field, bottom_plane)\n",
    "r = r1 + r2 + r3\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20.2.4 - Surface Integrals of Vector Fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rs is the symbol representing the radius\n",
    "rs, phi, theta = symbols(\"r_s, phi, theta\")\n",
    "r = 3\n",
    "semi_sphere = ParametricRegion(\n",
    "    (r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)),\n",
    "    (phi, 0, pi/2), (theta, 0, 2*pi)\n",
    ")\n",
    "disk = ParametricRegion(\n",
    "    (rs*cos(theta), rs*sin(theta), 0),\n",
    "    (theta, 0, 2*pi), (rs, 0, r)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from spb import *\n",
    "(s_x, s_y, s_z), s_r1, s_r2 = semi_sphere.args\n",
    "(d_x, d_y, d_z), d_r1, d_r2 = disk.args\n",
    "graphics(\n",
    "    # semi sphere\n",
    "    surface_parametric(s_x, s_y, s_z, s_r1, s_r2, wireframe=True, wf_n2=20),\n",
    "    # disk\n",
    "    surface_parametric(d_x, d_y, d_z, d_r1, d_r2, wireframe=True, wf_n1=20),\n",
    "    backend=PB, aspect=\"equal\", legend=False, grid=False\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "semi_sphere_negative = ParametricRegion(\n",
    "    (r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)),\n",
    "    (theta, 0, 2*pi), (phi, 0, pi/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = x*i + y*j + z**4 * k\n",
    "r1 = vector_integrate(F, semi_sphere)\n",
    "r2 = vector_integrate(F, disk)\n",
    "res = r1 + r2\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 20.2.5 - Volume Integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = CoordSys3D(\"C\")\n",
    "x, y, z = C.base_scalars()\n",
    "r, phi, theta = symbols(\"r, phi, theta\")\n",
    "volume = ParametricRegion(\n",
    "    (r*sin(phi)*cos(theta), r*sin(phi)*sin(theta), r*cos(phi)),\n",
    "    (r, 0, 2), (phi, 2*pi/3, pi), (theta, pi / 2, 3*pi/2)\n",
    ")\n",
    "vector_integrate(x * z, volume)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from spb import *\n",
    "\n",
    "var(\"r, phi\")\n",
    "phi_val = 2*pi/3\n",
    "r_val = 2\n",
    "\n",
    "graphics(\n",
    "    # cone\n",
    "    surface_parametric(\n",
    "        r*sin(phi_val)*cos(theta), r*sin(phi_val)*sin(theta), r*cos(phi_val),\n",
    "        (r, 0, r_val), (theta, pi / 2, 3*pi/2), wireframe=True\n",
    "    ),\n",
    "    # sphere cap\n",
    "    surface_parametric(\n",
    "        r_val*sin(phi)*cos(theta), r_val*sin(phi)*sin(theta), r_val*cos(phi),\n",
    "        (phi, 2*pi/3, pi), (theta, pi / 2, 3*pi/2), wireframe=True\n",
    "    ),\n",
    "    # sphere full\n",
    "    surface_parametric(\n",
    "        r_val*sin(phi)*cos(theta), r_val*sin(phi)*sin(theta), r_val*cos(phi),\n",
    "        (phi, 0, 2*pi/3), (theta, pi / 2, 3*pi/2),\n",
    "        wireframe=True, wf_n2=10, wf_n1=15,\n",
    "        wf_rendering_kw={\"line_color\": \"gray\", \"line_dash\": \"dot\", \"line_width\": 2},\n",
    "        rendering_kw={\"opacity\": 0.0}\n",
    "    ),\n",
    "    backend=PB, grid=False, legend=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region1 = ParametricRegion(\n",
    "    (x, y, z), (z, 0, 6 - 2*x - 3*y), (y, 0, -2*x/3 + 2), (x, 0, 3))\n",
    "vector_integrate(2*x, region1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region2 = ParametricRegion(\n",
    "    (x, y, z), (z, 0, 6 - 2*x - 3*y), (x, 0, -3*y/2 + 3), (y, 0, 2))\n",
    "vector_integrate(2*x, region2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate(2*x, (z, 0, 6 - 2*x - 3*y), (y, 0, -2*x/3 + 2), (x, 0, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from spb import *\n",
    "\n",
    "def plane_from_eq(expr, x, y, z):\n",
    "    z_expr = solve(expr, z)[0]\n",
    "    p1 = Point3D(0, 1, z_expr.subs({x: 0, y: 1}))\n",
    "    p2 = Point3D(1, 0, z_expr.subs({x: 1, y: 0}))\n",
    "    p3 = Point3D(1, 1, z_expr.subs({x: 1, y: 1}))\n",
    "    return Plane(p1, p2, p3)\n",
    "\n",
    "expr = 2*x + 3*y + z - 6\n",
    "plane = plane_from_eq(expr, x, y, z)\n",
    "p = plot_geometry(\n",
    "    plane, (x, 0, 5), (y, 0, 5), (z, 0, 7),\n",
    "    backend=PB, n=500, xlabel=\"x\", ylabel=\"y\", zlabel=\"z\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 20.2 - Advanced Topics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Click here to open the associated notebook](Chapter-20-Advanced-Topics.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}