#!/usr/bin/env python
# coding: utf-8

# ## 20. Integration of differential forms on singular $k$-cubes
# 
# This notebook is part of the [Introduction to manifolds in SageMath](https://sagemanifolds.obspm.fr/intro_to_manifolds.html) by Andrzej Chrzeszczyk (Jan Kochanowski University of Kielce, Poland).

# In[1]:


version()


# **Remark**. This notebook is an attempt to present a more friendly version of Chapter 4 of:
# 
# Michael Spivak. Calculus on manifolds, Addison Wesley, NY 1965,
# 
# with elementary examples in SageMath.
# 
# **Remark**. Readers interested in application of free software to illustrate  the classical Stokes theorem should consult for example https://docs.sympy.org/latest/modules/vector/vector_integration.html
# <br>
# 
# **Warning**. In the present notebook we extensively use the information on differential forms from [notebook 14](https://nbviewer.org/github/sagemanifolds/IntroToManifolds/blob/main/14Manifold_Differential_Forms.ipynb) and 
# [notebook 16](https://nbviewer.org/github/sagemanifolds/IntroToManifolds/blob/main/16Manifold_Ext_Der.ipynb).
# 
# <br>
# 
# ### Explicit description of some pullbacks
# 
# Let us start from some more explicit description of pullback of $k$-forms on $k$-dimensional open sets.
# 
# Consider the smooth map $\phi:U\to V$ between open subsets in $R^k$. Let $x^1,\ldots,x^k$ be the coordinates on $V$ and $t^1,\ldots,t^k$ - the coordinates on $U$. For $k$-form $\omega(x)=f(x)dx^1\wedge\ldots\wedge dx^k,\ $ $f\in C^\infty(V)$ 
# we have
# 
# \begin{equation}
# (\phi^*\omega)(t)=f(\phi(t))\frac{\partial(\phi^1,\ldots,\phi^k)}{\partial(t^{1},\ldots,t^{k})}dt^1\wedge\ldots\wedge dt^k,
# \tag{20.1}
# \end{equation}
# where $\phi=(\phi^1,\ldots,\phi^k)$ and $\ \frac{\partial(\phi^1,\ldots,\phi^k)}{\partial(t^{1},\ldots,t^{k})}=\det\big[\frac{\partial\phi^i}{\partial t^j}\big]$.
# 
# In fact
# 
# $$\phi^*(\omega)=f(\phi(t))\big(\frac{\partial\phi^1}{\partial t^{i_1}}dt^{i_1}\big)
# \wedge\ldots\wedge\big(\frac{\partial\phi^k}{\partial t^{i_k}}dt^{i_k}\big)\\
# =f(\phi(t))\frac{\partial\phi^1}{\partial t^{i_1}}\ldots\frac{\partial\phi^k}{\partial t^{i_k}}dt^{i_1}\wedge \ldots\wedge dt^{i_k}\\
# =f(\phi(t))\big(\sum_{\sigma\in S_k} (\text{sign}\,\sigma)
# \frac{\partial\phi^1}{\partial t^{\sigma(1)}}\ldots\frac{\partial\phi^k}{\partial t^{\sigma(k)}}\big)dt^{1}\wedge \ldots\wedge dt^{k}\\
# =f(\phi(t))\det\big[\frac{\partial\phi^i}{\partial t^j}\big]dt^{1}\wedge \ldots\wedge dt^{k},
# $$
# 
# (we have used the relations $i_1=\sigma(1),\ldots,i_k=\sigma(k),\ $ $\ dt^{i_1}\wedge\ldots\wedge dt^{i_k}=(\text{sign}\,\sigma)dt^1\wedge\ldots\wedge dt^k$ and the definition of determinant).

# <br>
# 
# ### Singular $k$-cubes
# 
# <br>
# 
# Let
# $$[0,1]^k=\underbrace{[0,1]\times\ldots\times[0,1]}_{k\mbox{-times}}$$
# and $U$ be an open subset of $R^k$ such that $[0,1]^k\subset U.$
# 
# Let $\phi : U \to V$
# be a $C^\infty$ map  of $U$ into an open set $V$
# in $R^m$ ($k$ and $m$ may be different). We will write this map simply as
# \begin{equation}
# \phi:[0,1]^k\to V,
# \tag{20.2}
# \end{equation}
# 
# (this map is a restriction of a  $C^\infty$ map $ U\to V$).  
# 
# The map $\phi:[0,1]^k\to V$ will be
# called a **singular $k$-cube** in $V$. The word singular indicates that the map need not be one-to-one,  though we will require that $\phi$ be in $C^\infty$. For $k=1,2$ the singular $k$-cubes are just parametrized curves and surfaces, respectively.
# 
# <br>
# 
# ### Standard $k$-cubes
# 
# <br>
# 
# By the **standard $k$-cube** we will mean the map
# $ I:[0,1]^k\to R^k$ defined by
# 
# \begin{equation}
# I(t)=t,\quad t\in[0,1]^k,
# \tag{20.3}
# \end{equation}
# 
# so $I$ is an example of singular $k$-cube with $\phi\big|_{[0.1]^k}=id_{[0,1]^k}.$
# 
# For standard $k$-cube $I$ we can define singular $(k-1)$-cubes $I^0_i$ and $I^1_i$
# (the i-th pieces of the boundary of $I$) as the maps
# 
# \begin{equation}
# \begin{matrix}
# I^0_i:[0,1]^{k-1}\ni t\to x(t)\in R^k, \text{ where}\\
# x^1(t)=t^1,\\
# ...........,\\
# x^{i-1}(t)=t^{i-1},\\
# x^i(t)=x^i=0,\\
# x^{i+1}(t)=t^i,\\
# ...........,\\
# x^k(t)=t^{k-1},\\
# t=(t^1,\ldots,t^{k-1})\in [0,1]^{k-1},
# \end{matrix}
# \tag{20.4}
# \end{equation}
# 
# and
# 
# \begin{equation}
# \begin{matrix}
# I^1_i:[0,1]^{k-1}\ni t\to x(t)\in R^k, \text{ where}\\
# x^1(t)=t^1,\\
# ...........,\\
# x^{i-1}(t)=t^{i-1},\\
# x^i(t)=x^i=1,\\
# x^{i+1}(t)=t^i,\\
# ...........,\\
# x^k(t)=t^{k-1},\\
# t=(t^1,\ldots,t^{k-1})\in [0,1]^{k-1}.
# \end{matrix}
# \tag{20.5}
# \end{equation}
# 
# 
# The singular $(k-1)$-cubes $I^0_i$ and  $I^1_i$ are called $(0,i)$ and $(1,i)$ faces of $I$ respectively.

# <br>
# 
# ### Boundary of the singular $k$-cube
# 
# <br>
# 
# In the definition of the boundary of a singular $k$-cube we will use some auxiliary notions.
# 
# For arbitrary family of objects $W=\{w_i\}_{i\in J}$ we can construct a **vector space of formal linear combinations of elements $w_i\in W$** as  the family of functions $f:W\to R$, such that $f(w_i) = 0$
# for all but finitely many $w_i$. Vector space operations on such functions are defined as usual: if $f,g$ are formal linear combinations of elements of $W$, then
# 
# $$(f+g)(w_i)=f(w_i)+g(w_i)\quad\mbox{and }\ \ (qf)(w_i)=q(f(w_i)),\quad \mbox{for } w_i\in W, q\in R.$$
# 
# The formal linear combination $f$ are usually written as 
# 
# $$
# f=\sum a_iw_i,\quad \mbox{ where } \ \ a_i=f(w_i).
# $$
# 
# If we take  as $W$ the family of all singular $k$-cubes in $V$ we can define the **space of $k$-chains in $V$** as the space of formal linear combinations  of singular $k$-cubes in $V$. Thus  the space of $k$-chains consists of (finite) linear combinations 
# 
# $$ \sum a_i\phi_i,\quad \mbox{where } \phi_i\ \ \mbox{are singular}\ \  \text{$k$-cubes and}\ \ a_i\in R.
# $$

# Using the introduced notion we can define the boundary $\partial I$ of the **standard $k$-cube $I$** as the $(k-1)$-chain
# 
# \begin{equation}
# \partial I=\sum_{i=1}^k (-1)^{i-1}(I^1_i-I^0_i).
# \tag{20.6}
# \end{equation}
# 
# In the case of **singular $k$-cube $\phi:[0,1]^k\to V$** we define the 
# $(0,i)$ and $(1,i)$ faces as
# 
# \begin{equation}
# \phi^0_i=\phi\circ I^0_i,\qquad \phi^1_i=\phi\circ I^1_i,
# \tag{20.7}
# \end{equation}
# 
# respectively and the boundary $\partial \phi$ of $\phi$ as the $(k-1)$-chain
# 
# \begin{equation}
# \partial\phi =\sum_{i=1}^k (-1)^{i-1}(\phi^1_i-\phi^0_i).
# \tag{20.8}
# \end{equation}

# ### Integrals of $k$-forms on singular $k$-cubes and chains
# 
# <br>

# Let $\omega$ be a smooth $k$-form on an open set $U\subset R^k$ containing $[0,1]^k$.
# Assume that $\omega$ is  of the form
# 
# \begin{equation}
# \omega=fdx^1\wedge\ldots\wedge dx^k,
# \tag{20.9}
# \end{equation}
# 
# where $f$ is a smooth function on $U$.
# 
# We define first the integral of $\omega$ over  $[0,1]^k$
# 
# \begin{equation}
# \int_{[0,1]^k}\omega =\int_{[0,1]^k}fdx^1\wedge\ldots\wedge dx^k=\int_{[0,1]^k}f(x)dx^1\ldots dx^k=\int_{[0,1]^k} f,
# \tag{20.10}
# \end{equation}
# 
# where the right hand side denotes the Riemann integral over the rectangle $[0,1]^k$.
# 
# Using this relation we can define  the **integral of the smooth $k$-form $\omega$ on $V$ over singular $k$-cube** $\phi:[0,1]^k\to V$:
# 
# \begin{equation}
# \int_\phi\omega =\int_{[0,1]^k}\phi^\star\omega.
# \tag{20.11}
# \end{equation}
# 
# With the help of (20.1) we can rewrite this formula: 
# 
# \begin{equation}
# \int_{\phi}\omega=\int_{[0,1]^k}f(\phi(t))
# \frac{\partial(\phi^1,\ldots,\phi^k)}{\partial(t^{1},\ldots,t^{k})}dt^1\wedge\ldots\wedge dt^k,
# \tag{20.12}
# \end{equation}
# 
# where $\phi=(\phi^1,\ldots,\phi^k)$.
# 
# **In the special case of standard $k$-cube $I$** given by $x^1=t^1,\ldots,x^k=t^k$, $t\in [0,1]^k$ we obtain
# 
# \begin{equation}
# \begin{matrix}
# \int_I fdx^1\wedge\ldots\wedge dx^k=
# \int_{[0,1]^k}I^*(fdx^1\wedge\ldots\wedge dx^k)\\
# =\int_{[0,1]^k}f(I(t))
# \frac{\partial(x^1,\ldots,x^k)}{\partial(t^{1},\ldots,t^{k})}dt^1\wedge\ldots\wedge dt^k\\
# =\int_{[0,1]^k}f(t)dt^1\ldots dt^k.
# \end{matrix}
# \tag{20.13}
# \end{equation}
# 
# Finally we define the **integral of $k$-form of $\omega$ over $k$-chain**  $c=\sum a_i\phi_i$:
# 
# \begin{equation}
# \int_c\omega=\sum a_i\int_{\phi_i}\omega.
# \tag{20.14}
# \end{equation}
# 
# The last formula applied to $(k-1)$-chain $\ \ \partial I=\sum_{i=1}^k (-1)^{i-1}(I^1_i-I^0_i)\ \ $ and $(k-1)$-form $\omega$ gives
# 
# \begin{equation}
# \int_{\partial I}\omega=\sum_{i=1}^k (-1)^{i-1}(\int_{I^1_i}\omega
# -\int_{I^0_i}\omega),
# \tag{20.15}
# \end{equation}
# 
# and in the case of $(k-1)$-chain $\ \ \partial \phi=\sum_{i=1}^k (-1)^{i-1}(\phi^1_i-\phi^0_i)\ \ $: 
# 
# \begin{equation}
# \int_{\partial \phi}\omega=\sum_{i=1}^k (-1)^{i-1}(\int_{\phi^1_i}\omega
# -\int_{\phi^0_i}\omega).
# \tag{20.16}
# \end{equation}

# <br>
# 
# 
# The definition (20.15) of the integral of $\omega$  on $\partial I$ looks artificially, but we can give it a quite natural interpretation if we define the orientation of $I^0_i$ and $I^1_i$.

# We recall  that the transition from one frame $e_1,..., e_k$ in $R^k$  to a
# second frame $e_1',..., e_k'$ can be defined by the square matrix obtained from  $e_j'=A^i_je_i$. The determinant of this matrix is always nonzero, and the set of all frames divides into two equivalence classes, each class containing all possible frames such that for any two of them the  determinant of  transition matrix is positive. Such equivalence classes are called orientation classes of frames in $R^n$.
# To define an orientation means to fix one of these orientation classes. Thus,  **the oriented space $R^k$ is the space $R^k$  together with a fixed orientation class of frames**.

# Let $n_i$ be the unit  outward (with respect to $[0,1]^k$) normal vector to $I^1_i$ or $I^0_i$ defined in (20.4),(20,5).
# Assume that the space $R^𝑘$ has been oriented by $e_1,\ldots,e_k$.
# We can define the **orientation of $I^1_0$ and $I^0_i$ prescribed by  the unit outward normal vector to $I^1_i$ or $I^0_i$** in the following way.
# 
# Compare the frame $𝑛_i,𝑒_1,...,\widehat{e_i}\ldots,𝑒_k$ ($e_i$ omitted) with the frame $e_1,..,,𝑒_k$, that orients the space $R^𝑘$. 
# 
# If these are in the same orientation class, then $I^0_i$ (or $I^1_i$)
# has the orientation prescribed by  the unit outward normal vector to $I^0_i$ (or $I^1_i$). 
# 
# If these two frames are in different orientation classes, then we take the orientation of $I^0_i$ ( or $I^1_i$)  opposite to that defined by $𝑛_i,𝑒_1,...,\widehat{e_i}\ldots,𝑒_k$ as the one prescribed by the outward  unit normal vector.
# 
# Defining $\ \int_{\partial I}\omega\ \ $ we take the sum of all $\ \int_{ I^0_i}\omega\ \ $ and 
# $\ \int_{I^1_i}\omega,\ \ $ but if $\ I^0_i\ $ or $\ I^1_i\ $ have the orientation opposite to that prescribed by the outward unit normal vector we take the corresponding integrals with minus sign. 

# <br>
# 
# **Example 20.1**
# 
# Consider the case $k=2$.

# In[2]:


l1=polygon2d([(1.0,0.5),(1.5,0.5),(1.5,1.0),(1.0,1.0)],
             fill=False,color='black',xmin=0,xmax=2.5,
            ymin=0,ymax=1.5)
e1=arrow((1.35,1),(1.1,1),color='black')
t1=text('$-e_1$',(1.1,1.1),fontsize=20, color='black')
n1=arrow((1.35,1),(1.35,1.2),color='black')
tn1=text('$n$',(1.35,1.25),fontsize=20, color='black')
I21=text('$I^1_2$',(1.2,1.3),fontsize=24, color='black')

e2=arrow((1.1,0.5),(1.35,0.5),color='black')
t2=text('$e_1$',(1.25,0.6),fontsize=20, color='black')
n2=arrow((1.1,0.5),(1.1,0.3),color='black')
tn2=text('$n$',(1.1,0.25),fontsize=20, color='black')
I20=text('$I^0_2$',(1.3,0.2),fontsize=24, color='black')

e3=arrow((1,0.85),(1.0,0.65),color='black')
t3=text('$-e_2$',(0.85,0.65),fontsize=20, color='black')
n3=arrow((1.0,0.85),(0.8,0.85),color='black')
nt3=text('$n$',(0.75,0.9),fontsize=20, color='black')
I11=text('$I^1_1$',(1.9,0.75),fontsize=24, color='black')

e4=arrow((1.5,0.65),(1.5,0.85),color='black')
t4=text('$e_2$',(1.6,0.85),fontsize=20, color='black')
n4=arrow((1.5,0.65),(1.7,0.65),color='black')
nt4=text('$n$',(1.75,0.65),fontsize=20, color='black')
I01=text('$I^0_1$',(0.6,0.75),fontsize=24, color='black')


(l1+e1+t1+n1+tn1+e2+n2+tn2+t2+e3+t3+n3+nt3+
 e4+t4+n4+nt4+I21+I20+I11+I01).show(axes=False)


# If $k=2$, then the integral  over the boundary of $\ I\ $  defined 
# in (20.15) takes the form
# 
# \begin{equation}
# \begin{matrix}
# \int_{\partial I}\omega
# =\sum_{i=1}^2(-1)^{i-1}\big[\int_{I^1_i}\omega-\int_{I^0_i}\omega\big] \\
# =\int_{I^1_1}\omega-\int_{I^0_1}\omega
# -\int_{I^1_2}\omega+\int_{I^0_2}\omega.
# \end{matrix}
# \tag{*}
# \end{equation}
# 
# The orientation of the boundary of $I$ prescribed by  the unit normal outward vector implies that: 
# 
# on $I^0_2$ the frame  $n,e_1$ is in the same orientation class as $e_1,e_2$ in $R^2$,
# 
# on $I^1_1$ the frame  $n,e_2$ is in the same orientation class as $e_1,e_2$ in $R^2$,
# 
# on $I^1_2$ the frame  $n,-e_1$ is in the same orientation class as $e_1,e_2$ in $R^2$,
# 
# on $I^0_1$ the frame  $n,-e_2$ is in the same orientation class as $e_1,e_2$ in $R^2$.
# 
# This is in accordance with the fact that in (*) the integrals over $I^1_2$ and $I^0_1$ are taken with the minus sign.
# 
# As we can see, in the case $k=2$ the formula (20.15) agrees with
# the orientation of boundary of $I$ prescribed by  the unit outward normal vector and  corresponds to going around the square
# counterclockwise.
# 

# <br>
# 
# **Example 20.2**
# 
# Consider the case $k=3$.

# In[3]:


x1=arrow((0.5,0.2),(1.8,0.2),color='black')
x2=arrow((0.5,0.2),(1.6,1.05),color='lightgrey')
x3=arrow((0.5,0.2),(0.5,1.4),color='black')
t1=text('$x_1$',(2.1,0.2),fontsize=20, color='black')
t2=text('$x_2$',(1.5,1.4),fontsize=20, color='black')
t3=text('$x_3$',(0.4,1.7),fontsize=20, color='black')
#(x1+x2+x3).show(axes=False)

l1=polygon2d([(0.5,0.2),(0.5,0.8),(1.0,0.8),(1,0.2)],
             fill=False,color='black')
l2=polygon2d([(0.5,0.8),(0.8,1.),(1.3,1.),(1.0,0.8)],
             fill=False,color='black')
l3=polygon2d([(1.0,0.2),(1.0,0.8),(1.3,1),(1.3,0.5)],
             fill=False,color='black')

p11=x1+t1+x2+t2+x3+t3+l1+l2+l3

e2=arrow((0,0),(-0.6,-0.6),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t2=text('$-e_2$',(-0.75,-0.75),fontsize=20, color='black')
e3=arrow((0,0),(0,1.2),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t3=text('$e_3$',(0,1.4),fontsize=20, color='black')
n=arrow((0,0),(-0.7,0),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
tn=text('$n$',(-0.8,0.),fontsize=20, color='black')
left=text("$I^0_1$",(1.8,0),fontsize=28, color='black')
p10=e2+t2+e3+t3+n+tn+left

e2=arrow((0,0),(0.6,0.6),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t2=text('$e_2$',(0.75,0.75),fontsize=20, color='black')
e3=arrow((0,0),(0,1.2),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t3=text('$e_3$',(0,1.4),fontsize=20, color='black')
n=arrow((0,0),(0.7,0),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
tn=text('$n$',(0.8,0.),fontsize=20, color='black')
right=text("$I^1_1$",(-1.9,0),fontsize=28, color='black')
p12=e2+t2+e3+t3+n+tn+right

e1=arrow((0,0),(0.7,0.0),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t1=text('$e_1$',(0.8,0.0),fontsize=20, color='black')
e2=arrow((0,0),(0.7,0.7),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t2=text('$e_2$',(0.8,0.8),fontsize=20, color='black')
n=arrow((0,0),(0,1),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
tn=text('$n$',(0.,1.1),fontsize=20, color='black')
ttop=text("$I^1_3$",(0,-1.9),fontsize=28, color='black')
p01=e1+t1+e2+t2+n+tn+ttop

e1=arrow((0,0),(0.7,0.0),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t1=text('$e_1$',(0.8,0.0),fontsize=20, color='black')
e2=arrow((0,0),(-0.7,-0.7),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t2=text('$-e_2$',(-0.8,-0.8),fontsize=20, color='black')
n=arrow((0,0),(0,-1),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
tn=text('$n$',(0.,-1.1),fontsize=20, color='black')
bottom=text("$I^0_3$",(0,1.8),fontsize=28, color='black')
p21=e1+t1+e2+t2+n+tn+bottom

e1=arrow((0,0),(-0.7,0.0),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t1=text('$-e_1$',(-1,0.0),fontsize=20, color='black')
e3=arrow((0,0),(0,1.2),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t3=text('$e_3$',(0,1.4),fontsize=20, color='black')
n=arrow((0,0),(0.7,0.7),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
tn=text('$n$',(0.8,0.8),fontsize=20, color='black')
back=text("$I^1_2$",(-1.8,-1.8),fontsize=28, color='black')
p02=e1+t1+e3+t3+n+tn+back

e1=arrow((0,0),(0.7,0.0),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t1=text('$e_1$',(0.8,0.0),fontsize=20, color='black')
e3=arrow((0,0),(0,1.2),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
t3=text('$e_3$',(0,1.4),fontsize=20, color='black')
n=arrow((0,0),(-0.7,-0.7),color='grey',xmin=-2,ymin=-2,xmax=2,ymax=2)
tn=text('$n$',(-0.8,-0.8),fontsize=20, color='black')
front=text("$I^0_2$",(1.8,1.8),fontsize=28, color='black')
p20=e1+t1+e3+t3+n+tn+front

graphics_array([[text("",(1,1)),p01,p02],
                [p10,p11,p12],
                [p20,p21,text("",(1,1))]]).show(axes=False,
                                    figsize=(8,6))


# If $k=3$, then the integral  over the boundary of $\ I\ $  defined 
# in (20.15) takes the form
# 
# \begin{equation}
# \begin{matrix}
# \int_{\partial I}\omega
# =\sum_{i=1}^3(-1)^{i-1}\big[\int_{I^1_i}\omega-\int_{I^0_i}\omega\big] \\
# =\int_{I^1_1}\omega-\int_{I^0_1}\omega
# -\int_{I^1_2}\omega+\int_{I^0_2}\omega
# +\int_{I^1_3}\omega-\int_{I^0_3}\omega.
# \end{matrix}
# \tag{**}
# \end{equation}
# 
# The orientation of the boundary of $I$ prescribed by  the unit outward normal vector implies that: 
# 
# on $I^1_1$ the frame  $n,e_2,e3$ is in the same orientation class as $e_1,e_2,e_3$ in $R^3$,
# 
# on $I^0_1$ the frame  $n,-e_2,e_3$ is in the same orientation class as $e_1,e_2,e_3$ in $R^3$,
# 
# on $I^1_2$ the frame  $n,-e_1,e_3$ is in the same orientation class as $e_1,e_2,e_3$ in $R^3$,
# 
# on $I^0_2$ the frame  $n,e_1,e_3$ is in the same orientation class as $e_1,e_2,e_3$ in $R^3$,
# 
# on $I^1_3$ the frame  $n,e_1,e_2$ is in the same orientation class as $e_1,e_2,e_3$ in $R^3$,
# 
# on $I^0_3$ the frame  $n,e_1,-e_2$ is in the same orientation class as $e_1,e_2,e_3$ in $R^3$.
# 
# This is in accordance with the fact that in (**) the integrals over $I^0_1,I^1_2,I^0_3$ are taken with the minus sign.
# 
# As we can see, in the case $k=3$ the formula (20.15) agrees with
# the orientation of boundary of $I$ prescribed by  the unit outward normal vector.
# 

# <br>
# 
# ### Integrals of $(k-1)$-forms over $I^0_i$ and $I^1_i$
# 
# <br>
# 
# 
# General  smooth $(k-1)$-form on an open subset $U\subset R^k$ containing  $[0,1]^k$ is of the form
# 
# \begin{equation}
# \omega=\sum_{j=1}^k\omega_j(x^1,\ldots,x^k)dx^1\wedge\ldots \widehat{dx^j}
# \wedge\ldots dx^k,
# \tag{20.17}
# \end{equation}
# ($dx^j$ is omitted). 
# 
# Consider the  single summand $\ \omega_j
# dx^1\wedge\ldots \wedge \widehat{dx^j}\wedge\ldots\wedge dx^k$.
# 
# Since the map  $I^0_i$ is
# given by (20.4)
# and $x^i=0$ is constant, for $i\neq j$ we have
# 
# $$
# \int_{I^0_i}\omega_j
# dx^1\wedge\ldots \wedge \widehat{dx^j}\wedge\ldots\wedge dx^k\\
# =\int_{[0,1]^{k-1}}(I^0_i)^*(\omega_j
# dx^1\wedge\ldots \wedge \widehat{dx^j}\wedge\ldots\wedge dx^k)\\
# =\int_{[0,1]^{k-1}}\omega_j(I^0_i(t))
# \frac{\partial(x^1,\ldots 0\ldots\widehat{x^j}\ldots,x^k)}{\partial(t^{1},\ldots,t^{k-1})}
# dt^{1}\wedge\ldots\wedge dt^{k-1}=0
# $$
# 
# and analogously for $I^1_i$ , so if $i\neq j$
# 
# \begin{equation}
# \begin{matrix}
# \int_{I^0_i}\omega_j
# dx^1\wedge\ldots \wedge \widehat{dx^j}\wedge\ldots\wedge dx^k=0,\\
# \int_{I^1_i}\omega_j
# dx^1\wedge\ldots \wedge \widehat{dx^j}\wedge\ldots\wedge dx^k=0.
# \end{matrix}
# \tag{20.18}
# \end{equation}
# 
# If $i=j$ we have
# 
# $$\int_{I^0_i}\omega_i
# dx^1\wedge\ldots \wedge \widehat{dx^i}\wedge\ldots\wedge dx^k\\
# =\int_{[0,1]^{k-1}}(I^0_i)^*(\omega_i
# dx^1\wedge\ldots \wedge \widehat{dx^i}\wedge\ldots\wedge dx^k)\\
# =\int_{[0,1]^{k-1}}\omega_i(I^0_i(t))
# \frac{\partial(x^1,\ldots\widehat{x^i}\ldots,x^k)}{\partial(t^{1},\ldots,t^{k-1})}
# dt^{1}\wedge\ldots\wedge dt^{k-1}\\
# =\int_{[0,1]^{k-1}}\omega_i(I^0_i(t))
# dt^{1}\ldots dt^{k-1}\\
# =\int_{[0,1]^{k-1}}\omega_i\big|_{x^i=0}
# dx^{1}\ldots\widehat{dx^i}\ldots dx^{k}\\
# =\int_{[0,1]^{k}}\omega_i\big|_{x^i=0}
# dx^{1}\ldots\ldots dx^{k}.
# $$
# 
# The last equality follows from the fact, that if a function $g$ is constant with respect to $x^i$, then $g=\int_0^1 gdx^i$.
# 
# 
# Analogously we can compute the integral over $I^1_i$, so
# 
# \begin{equation}
# \begin{matrix}
# \int_{I^0_i}\omega_i
# dx^1\wedge\ldots \wedge \widehat{dx^i}\wedge\ldots\wedge dx^k\\
# =\int_{[0,1]^{k}}\omega_i\big|_{x^i=0}
# dx^{1}\ldots dx^{k},
# \end{matrix}
# \tag{20.19}
# \end{equation}
# and
# \begin{equation}
# \begin{matrix}
# \int_{I^1_i}\omega_i
# dx^1\wedge\ldots \wedge \widehat{dx^i}\wedge\ldots\wedge dx^k\\
# =\int_{[0,1]^{k}}\omega_i\big|_{x^i=1}
# dx^{1}\ldots dx^{k}.
# \end{matrix}
# \tag{20.20}
# \end{equation}

# <br>
# 
# ### Stokes theorem on the standard $k$-cube
# 
# <br>
# 
# Consider a smooth $(k-1)$-form (20.17) on an open subset $V\subset R^k$ containing $[0,1]^k$ and
#  let $I$ be the standard $k$-cube defined in (20.3).
# Let us compute $\ \int_I d\omega.$
# 
# By linearity we can compute each summand
# 
# $$\int_I d(\omega_idx^1\wedge\ldots \widehat{dx^i}\wedge\ldots dx^k)$$
# 
# separately. Let us focus on $i=1$ first. We have
# 
# $$\int_I d(\omega_1 \widehat{dx^1}\wedge dx^2\wedge\ldots\wedge dx^k)\\
# =\int_{[0,1]^k}I^*( d(\omega_1 dx^2\wedge\ldots\wedge dx^k))\\
# =\int_{[0,1]^k} d(I^*(\omega_1 dx^2\wedge\ldots\wedge dx^k))\\
# =\int_{[0,1]^k}\frac{\partial \omega_1}{\partial t^1}dt^1\wedge\ldots\wedge dt^k\\
# =
# \int_{t^k=0}^{1}\ldots\int_{t^1=0}^{1}
# \frac{\partial \omega_1}{\partial t^1}dt^1\ldots dt^k\\
# =\int_{t^k=0}^{1}\ldots\int_{t^2=0}^{1}
# \big[\omega_1(t)\big|_{t^1=1}-\omega_1(t)\big|_{t^1=0}\big]dt^2\ldots dt^k. 
# $$
# 
# Since the integrand $g=\omega_1(t)\big|_{t^1=1}-\omega_1(t)\big|_{t^1=0}$ does not depend on $t^1$ it is equal to its integral $\int^1_0g\;dt^1$, so
# 
# $$\int_I d(\omega_1 \widehat{dx^1}\wedge dx^2\wedge\ldots\wedge dx^k)\\
# =\int_{[0,1]^k}\big[\omega_1(t)\big|_{t^1=1}-\omega_1(t)\big|_{t^1=0}\big]dt^1\ldots dt^k\\
# =\int_{[0,1]^k}\big[\omega_1\big|_{x^1=1}-\omega_1\big|_{x^1=0}\big]dx^1\ldots dx^k\\
# =\int_{I_1^1}\omega_1\widehat{dx^1}\wedge dx^2\wedge\ldots\wedge dx^k
# -\int_{I_1^0}\omega_1\widehat{dx^1}\wedge dx^2\wedge\ldots\wedge dx^k.
# $$
# 
# In the last equality we have used (20.19) and (20.20) for $i=1$.

# 
# 
# For $i=2$ we have
# 
# $$\int_I d(\omega_2 dx^1\wedge\widehat{dx^2}\wedge dx^3\ldots\wedge dx^k)\\
# =\int_{[0,1]^k}I^*( d(\omega_2 dx^1\wedge dx^3\ldots\wedge dx^k))\\
# =\int_{[0,1]^k} d(I^*(\omega_2 dx^1\wedge dx^3\ldots\wedge dx^k))\\
# =\int_{[0,1]^k}\frac{\partial \omega_2}{\partial t^2}dt^2\wedge dt^1\wedge dt^3\ldots\wedge dt^k\\
# =(-1)\int_{[0,1]^k}\frac{\partial \omega_2}{\partial t^2}dt^1\wedge dt^2\wedge dt^3\ldots\wedge dt^k\\
# =(-1)\int_{t^k=0}^{1}\ldots\int_{t^2=0}^1\int_{t^1=0}^{1}
# \frac{\partial \omega_2}{\partial t^2}dt^1\ldots dt^k\\
# =(-1)\int_{t^k=0}^{1}\ldots\int_{t^3=0}^1\int_{t^1=0}^1\int_{t^2=0}^{1}
# \frac{\partial \omega_2}{\partial t^2}dt^2dt^1 dt^3\ldots dt^k\\
# =(-1)\int_{t^k=0}^{1}\ldots\widehat{\int_{t^2=0}^1}\int_{t^1=0}^{1}
# \big[\omega_2(t)\big|_{t^2=1}-\omega_2(t)\big|_{t^2=0}\big]dt^1\ldots\widehat{dt^2}\ldots dt^k. 
# $$
# 
# Since the integrand $g=\omega_2(t)\big|_{t^2=1}-\omega_2(t)\big|_{t^2=0}$ does not depend on $t^2$ it is equal to its integral $\int^1_0g\;dt^2$, so
# 
# $$\int_I d(\omega_2 dx^1\wedge dx^3\wedge\ldots\wedge dx^k)\\
# =(-1)\int_{[0,1]^k}\big[\omega_2(t)\big|_{t^2=1}-\omega_2(t)\big|_{t^2=0}\big]dt^1\ldots dt^k\\
# =(-1)\int_{[0,1]^k}\big[\omega_2\big|_{x^2=1}-\omega_2\big|_{x^2=0}\big]dx^1\ldots dx^k\\
# =(-1)\big[\int_{I_2^1}\omega_2dx^1\wedge\widehat{dx^2}\wedge\ldots\wedge dx^k
# -\int_{I_2^0}\omega_2dx^1\wedge\widehat{dx^2}\wedge\ldots\wedge dx^k\big].
# $$
# 
# In the last equality we have used (20.19) and (20.20) for $i=2$.

# Now we are ready to consider the complete $(k-1)$ form (20.17). First compute the exterior differential
# 
# $$d\omega=\sum_{i=1}^k\frac{\partial \omega_i}{\partial x^i}
# dx^i\wedge dx^1\wedge\ldots \widehat{dx^i}\wedge\ldots\wedge dx^k\\
# =\sum_{i=1}^k (-1)^{i-1}\frac{\partial \omega_i}{\partial x^i}
# dx^1\wedge\ldots \wedge dx^i\wedge\ldots\wedge dx^k\\
# =\big(\sum_{i=1}^k (-1)^{i-1}\frac{\partial \omega_i}{\partial x^i}\big)
# dx^1\wedge\ldots \wedge dx^i\wedge\ldots\wedge dx^k,
# $$
# 
# and next the integral
# 

# $$\int_{I}d\omega=\sum_{i=1}^k(-1)^{i-1}\int_{I}
# \frac{\partial \omega_i}{\partial x^i}
# dx^1\wedge\ldots \wedge dx^i\wedge\ldots\wedge dx^k\\
# =\sum_{i=1}^k(-1)^{i-1}\int_{t^k=0}^{1}
# \ldots\widehat{\int_{t^i=0}^{1}}\ldots\int_{t^1=0}^{1}\int_{t^i=0}^{1}
# \frac{\partial \omega_i}{\partial t^i}
# dt^idt^1\ldots \widehat{dt^i}\ldots dt^k\\
# =\sum_{i=1}^k(-1)^{i-1}\int_{t^k=0}^{1}
# \ldots\widehat{\int_{t^i=0}^{1}}\ldots\int_{t^1=0}^{1}
# [\omega_i\big|_{x^i=1}-\omega_i\big|_{x^i=0}]dt^1\ldots \widehat{dt^i}\ldots dt^k\\
# =\sum_{i=1}^k(-1)^{i-1}\int_{[0,1]^k}
# \big[\omega_i\big|_{x^i=1}-\omega_i\big|_{x^i=0}\big]dx^1\ldots \ldots dx^k\\
# =\sum_{i=1}^k(-1)^{i-1}\big[\int_{I^1_i}\omega_i
# dx^1\wedge\ldots \wedge \widehat{dx^i}\wedge\ldots\wedge dx^k\\
# -\int_{I^0_i}\omega_i
# dx^1\wedge\ldots \wedge \widehat{dx^i}\wedge\ldots\wedge dx^k\big].
# $$

# As a consequence of (20.18)-(20.20) we obtain
# 
# $$
# \int_{I^0_i}\omega_i
# dx^1\wedge\ldots \wedge \widehat{dx^i}\wedge\ldots\wedge dx^k
# =\int_{I^0_i}\omega,\\
# \int_{I^1_i}\omega_i
# dx^1\wedge\ldots \wedge \widehat{dx^i}\wedge\ldots\wedge dx^k
# =\int_{I^1_i}\omega.
# $$
# 
# Thus we have checked that
# 
# \begin{equation}
# \int_{I}d\omega=
# \sum_{i=1}^k(-1)^{i-1}\big[\int_{I^1_i}\omega-\int_{I^0_i}\omega\big].
# \tag{20.21}
# \end{equation}

# Recall, that in (20.15) we defined the boundary integral of $\omega$ as
# $
# \int_{\partial I}\omega=\sum_{i=1}^k(-1)^{i-1}\big[\int_{I^1_i}\omega-\int_{I^0_i}\omega\big],
# $
# so we get:
# 
# <br>
# 
# ### Formulation of  Stokes theorem on the standard $k$-cube.
# 
# <br>
# 
# If $\ \omega\ $ is a smooth $(k-1)$-form (20.17) on an open set $U\subset R^k$ containing $[0,1]^k$,$\ \ I$ is the standard $k$-cube defined in (20.3)
# and if the boundary integral $\ \int_{\partial I}\omega\ $ is defined by  (20.15), then
# 
# \begin{equation}
# \int_{I}d\omega=\int_{\partial I}\omega.
# \tag{20.22}
# \end{equation}

# <br>
# 
# **Example 20.3**
# 
# Using Stokes theorem compute the integral
# $\ \int_{\partial I} xy^2zdx\wedge dy,\ $ where $\partial I$ denotes the boundary of the unit cube
# $[0,1]\times[0,1]\times[0,1]$, oriented by the outward unit normal.
# 
# <br>
# 
# Using Stokes theorem we can replace the integral by $\int_I da,\ $ where $a=xy^2zdx\wedge dy,$
# 
# so we start from definition of 2-form $a$ and computing its exterior derivative $da$.

# In[4]:


get_ipython().run_line_magic('display', 'latex')
M = Manifold(3, 'R^3')            # manifold M=R^3
c_xyz.<x,y,z> = M.chart()         # Cartesian coordinates
a=M.diff_form(2,name='a')         # 2-form a
E=c_xyz.coframe()                 # coframe (dx,dy,dz)
dx,dy,dz=E[:]                     # define dx,dy,dz
a=x*y^2*z*dx.wedge(dy)            # define 2-form a
da=a.exterior_derivative()        # exterior derivative of a
da.disp()                         # show da


# Computing $\int_I da,\ $ we use (20.13) and the iterated integral
# $\int_0^1\int_0^1\int_0^1 xy^2dzdydx.$

# In[5]:


integrate(integrate(integrate(x*y^2,(z,0,1)),(y,0,1)),(x,0,1))


# <br>
# 
# ### Formulation  of Stokes theorem on singular $k$-cubes
# 
# <br>
# 
# Assume that $V$ is an open subset of $R^k$, let
# $\omega$ be a smooth $(k-1)$-form on $V$ and $\phi:I\to V$ 
# be singular $k$-cube.
# Recall that in (20.16) we defined the boundary integral of $\omega$ over $\partial \phi$ as
# $\ \ \int_{\partial\phi}\omega=
# \sum_{i=1}^k(-1)^{i-1}\big[\int_{\phi^1_i}\omega-\int_{\phi^0_i}\omega\big]$.
# 
# 
# The Stokes theorem on singular $k$-cubes reads as follows
# 
# \begin{equation}
# \int_{\phi} d\omega=
# \int_{\partial\phi}\omega.
# \tag{20.23}
# \end{equation}

# <br>
# 
# 
# Before the proof we  check that for arbitrary singular $k$-cube $\phi:I\to V$ and $(k-1)$-form on $V$ the  following relation holds true
# 
# \begin{equation}
# \int_{\partial\phi}\omega=\int _{\partial I}\phi^*\omega.
# \tag{20.24}
# \end{equation}
# 
# Let us compute both sides separately.
# 
# $$\int_{\partial\phi}\omega=\sum_{i=1}^k(-1)^{i-1}\big(\int_{\phi^1_i}\omega
# -\int_{\phi^0_i}\omega\big)\\
# =
# \sum_{i=1}^k(-1)^{i-1}\big(\int_{[0,1]^{k-1}}(\phi^1_i)^*\omega
# -\int_{[0,1]^{k-1}}(\phi^0_i)^*\omega\big).
# $$
# 
# From the properties of pullback we know, that $\phi^1_i=\phi\circ I^1_i$ implies $(\phi^1_i)^*=(I^1_i)^*\circ\phi^*$ and analogously $(\phi^0_i)^*=(I^0_i)^*\circ\phi^*$.    Thus
# 
# $$\int_{\partial\phi}\omega=
# \sum_{i=1}^k(-1)^{i-1}\big(\int_{[0,1]^{k-1}}(I^1_i)^*\circ\phi^*\omega
# -\int_{[0,1]^{k-1}}(I^0_i)^*\circ\phi^*\omega\big).
# $$
# 
# On the other hand from (20.15)
# 
# $$
# \int_{\partial I}\phi^*\omega=\sum_{i=1}^k (-1)^{i-1}\big(\int_{I^1_i}\phi^*\omega
# -\int_{I^0_i}\phi^*\omega\big)\\
# =\sum_{i=1}^k (-1)^{i-1}\big(\int_{[0,1]^{k-1}}(I^1_i)^*\circ\phi^*\omega
# -\int_{[0,1]^{k-1}}(I^0_i)^*\circ\phi^*\omega\big),
# $$
# 
# so both sides of (20.24) are equal.
# 
# 
# 
# The formula (20.23) is a consequence of (20.22) and (20.24), since
# 
# 
# $$\int_{\phi} d\omega=\int_I\phi^*d\omega
# =\int_Id(\phi^*\omega)=\int_{\partial I}\phi^*\omega=\int_{\partial\phi}\omega.
# $$
# 

# <br>
# 
# ### Extending the integration domains
# 
# <br>

# Assume that $f,g$ are smooth functions on $[a,b],$ $\ f\leq g\ $and define  the set
# 
# \begin{equation}
# D =
# \{(u, v) ∈ R^2 : a ≤ u ≤ b, f (u) ≤ v ≤ g(u)\}.
# \tag{20.25}
# \end{equation} 
# 
#  Consider the mapping
# 
# $$
# ψ: [0,1]^2\ni
# (s, t)  \to \big((1-s)a+sb,(1-t)f(1-s)a+sb)+tg((1-s)a+sb)\big)\in D\subset R^2.
# $$
# 
# One can prove that
# 
# $ ψ : [0,1]^2 → D$ is one-to-one and $\psi(I^2)=D$.
# 
# 
# 
# If $\phi:D\to R^2$ is a smooth map, then $φ ◦ ψ : [0,1]^2 → \phi(D)\subset R^2$ is a singular 2-cube, so for 1-forms on some neighborhood of $\phi(D)$  Stokes formula takes the form
# 
# $$\int_{\phi\circ\psi}d\omega=
# \int_{\partial(\phi\circ\psi)}\omega.$$
# 
# If we define $\ \ \int_{\phi_D}d\omega=\int_{\phi\circ\psi}d\omega\ \ $
# and $\ \ \int_{\partial \phi_D}\omega=\int_{\partial(\phi\circ\psi)}\omega,\ \ $ then
# 
# $$\int_{\phi_D}d\omega=\int_{\phi\circ\psi}d\omega=
# \int_{\partial(\phi\circ\psi)}\omega=\int_{\partial \phi_D}\omega,
# $$
# 
# so we can use Stokes theorem for domains of type (20.25).
# 
# 
# Thus Stokes theorem can be applied to domains  whose parameterizations can be reparameterized into unit cubes.
# 
# **Remark.** In practice, computing  integrals $\int_{\phi_D} \omega$ for $\omega$ of the form 
# $\omega=hdx^1\wedge dx^2$ on $\phi(D)$ we do not use reparametrizations but  the formula
# $$\int_{\phi_D}\omega=\int_a^b\int_{f(u)}^{g(u)}h(\phi(u,v))\frac{\partial(\phi^1,\phi^2)}{\partial(u,v)}dudv.$$

# One can prove that  Stokes theorem can be applied to analogous domains in higher dimensions, for example
# 
# $$D =
# \{(u_1,u_2,u_3 ) ∈ R^3
# : a_1 ≤ u_1 ≤ b_1, c(u_1)\leq u_2\leq d(u_1), f (u_1,u_2) ≤ u_3 ≤ g(u_1,u_2)\},
# $$
# 
# for some smooth functions $c,d,f,g$.

# <br>
# 
# **Example 20.4**
# 
# Compute $\ \int_{S^2}xdy\wedge dz+ydz\wedge dx+zdx\wedge dy,$
# where $S^2$ denotes the 2-dimensional unit sphere,  which can be considered as the (reparametrized version of) singular 2-cube $\phi:[0,\pi]\times[0,2\pi]\to R^3$
# 
# $$(\theta,\phi)\to (\sin(\theta)\cos(\phi), \sin(\theta)\sin(\phi),\cos(\theta)).$$
# 
# <br>
# 
# We define the 2-form $a$ first.

# In[6]:


get_ipython().run_line_magic('display', 'latex')
M = Manifold(3, 'R^3')            # manifold M=R^3
c_xyz.<x,y,z> = M.chart()         # Cartesian coordinates
N = Manifold(2, 'N')              # manifold N=S^2
c_sph.<theta,phi>=N.chart()       # spherical coordinates
# embedding S^2->R^3 defines the singular 2-cube phi
Phi = N.diff_map(M, (sin(theta)*cos(phi), sin(theta)*sin(phi),cos(theta)),
                 name='ph',latex_name=r'\Phi')
a=M.diff_form(2,name='a')         # 2-form on N  a
E=c_xyz.coframe()                 # coframe (dx,dy,dz)
dx,dy,dz=E[:]                     # define dx,dy,dz
a=x*dy.wedge(dz)+y*dz.wedge(dx)+z*dx.wedge(dy) # define a
a.disp()                          # show a


# To use (reparametrized version of) (20.11) we need the pullback $\Phi^*a$.

# In[7]:


Phi.pullback(a).disp()


# To obtain the final result we use the iterated integral
# $\int_0^{2\pi}\int_0^\pi \sin\theta\; d\theta d\phi$.

# In[8]:


integrate(integrate(sin(theta), (theta,0,pi)), (phi,0,2*pi))


# <br>
# 
# **Example 20.5**
# 
# Compute $\ \int_S \frac{1}{3}zdx\wedge dy\ $ where $S$ is the part of the elliptic cone 
# $\ \frac{x^2}{6}+\frac{y^2}{9}= z^2,\ $  bounded by the planes $z=-1,z=1$, which may be considered as (reparametrized version of) the singular 2-cube $\psi:[-1,1]\times[0,2\pi]\to R^3$
# 
# $$(Z,\phi)\to (2Z\cos(\phi), 3Z\sin(\phi),Z).
# $$
# 
# <br>
# 
# Let us define 2-form $a$ and its pullback under $\psi$.

# In[9]:


get_ipython().run_line_magic('display', 'latex')
M = Manifold(3, 'R^3')            # manifold M=R^3
c_xyz.<x,y,z> = M.chart()         # Cartesian coordinates
N = Manifold(2, 'N')              # manifold N=S^2
c_sph.<phi,Z>=N.chart()       # spher. coordinates
E=c_xyz.coframe()                 # coframe (dx,dy,dz)
dx,dy,dz=E[:]                     # define dx,dy,dz
a=1/3*z*dx.wedge(dy)              # 2-form a
# psi defines the singular 2-cube
psi = N.diff_map(M, (2*Z*cos(phi), 3*Z*sin(phi),Z),
                 name='psi',latex_name=r'\psi')
psi.pullback(a).disp()            # pullback psi^*a


# We use (reparametrized version of) (20.11) and replace the double integral by the iterated integral $\int_0^{2\pi}\int_{-1}^1(-2Z^2)\;dZd\phi.$

# In[10]:


integrate(integrate(-2*Z^2,(Z,-1,1)),(phi,0,2*pi))


# <br>
# 
# **Example 20.6**
# 
# Compute the integral $\int_S xdy\wedge dz +zdz\wedge dx$, where $S$ is the sphere
# $x^2+y^2+z^2=9.$
# 
# <br>
# 
# First we define the 2-form $a$,

# In[11]:


E.<x1,x2,x3>=EuclideanSpace(3,"E")           # Euclidean space E^3
cart=E.cartesian_coordinates()               # Cartesian coordinates
spher.<r,th,ph> = E.spherical_coordinates()  # spherical coordinates
Phi=E.coord_change(spher, cart)     # transition spher -> cart 
a=E.diff_form(2,name='a')           # 2-form a
cf=c_xyz.coframe()                  # coframe (dx,dy,dz)
dx,dy,dz=cf[:]                      # define dx,dy,dz
a=x*dy.wedge(dz)+z*dz.wedge(dx)     # define 2-form a
a.disp()                            # show a


# and its exterior differential:

# In[12]:


da=a.exterior_derivative()
da.disp()


# Using spherical coordinates $(\theta,\phi)$, the sphere can be considered as reparametrized singular $2$-cube
# $\ \ [0,\pi]\times[0,2\pi]\to R^3$.

# In[13]:


Phi.disp()


# In[14]:


jac=Phi.jacobian_det();jac


# So if we apply (reparametrized version of) (20.12), and the relation $\int_0^{2\pi}d\phi=2\pi$, we obtain the following form of the iterated integral: 
# $2\pi(\int_0^3r^2\;dr)
# \int_0^\pi\sin(\theta)\;d\theta.$

# In[15]:


2*pi*integrate(r^2,(r,0,3))*integrate(sin(th),(th,0,pi))


# <br>
# 
# **Example 20.7**
# 
# Using Stokes theorem, compute the integral $\int_Sxdy\wedge dz+ydz\wedge dx + zdx\wedge dy$ over the unit sphere
# $x^2+y^2+z^2=1$. (In Example 20.4 we computed the integral using (20.11)).
# 
# <br>
# 
# As in the previous example we start from defining the 2-form $a$ and its exterior derivative.

# In[16]:


get_ipython().run_line_magic('display', 'latex')
M = Manifold(3, 'R^3')            # manifold M=R^3
c_xyz.<x,y,z> = M.chart()         # Cartesian coordinates
N = Manifold(2, 'N')              # manifold N=S^2
c_sph.<theta,phi>=N.chart()       # spherical coordinates
                                  # embedding S^2->R^3
psi = N.diff_map(M, (sin(theta)*cos(phi), sin(theta)*sin(phi),cos(theta)),
                 name='psi',latex_name=r'\psi')
a=M.diff_form(2,name='a')
E=c_xyz.coframe()                 # coframe (dx,dy,dz)
dx,dy,dz=E[:]                     # define dx,dy,dz
a=x*dy.wedge(dz)+y*dz.wedge(dx)+z*dx.wedge(dy)  # a
da=a.exterior_derivative()        # da
da.disp()


# We can use the value of Jacobian from previous example,

# In[17]:


# from previous example
jac=r^2*sin(theta)


# and obtain the following  iterated integral: 
# $\int_0^1\int_0^\pi\int^{2\pi}3\text{jac}\;d\phi d\theta dr.$

# In[18]:


integrate(integrate(integrate(3*jac,(phi,0,2*pi)),(theta,0,pi)),(r,0,1))


# <br>
# 
# **Example 20.8**
# 
# Compute the integral $ \int_{\partial S} xdx+xdy+ydz$, where $\partial S$ is the boundary of hemisphere $x^2+y^2+z^2=r^2,\quad z\geq 0$.
# 
# <br>
# 
# Define the 1-form $a$ and 2-form $da$.

# In[19]:


get_ipython().run_line_magic('display', 'latex')
M = Manifold(3, 'R^3')            # manifold M=R^3
c_xyz.<x,y,z> = M.chart()         # Cartesian coordinates
a=M.diff_form(2,name='a')
E=c_xyz.coframe()                 # coframe (dx,dy,dz)
dx,dy,dz=E[:]                     # define dx,dy,dz
a=x*dx+x*dy+y*dz                  # a
da=a.exterior_derivative()        # da
da.disp()


# Let us compute both sides of Stokes formula. Note that using spherical coordinates we can define the hemisphere as the (reparametrized) singular 2-cube
# 
# $$\Phi:[0,\frac{\pi}{2}]\times[0,2\pi]\ni(\theta,\phi)\to 
# (r\sin(\theta)\cos(\phi), r\sin(\theta)\sin(\phi),r\cos(\theta)).
# $$
# 
# The corresponding parts of $\partial\Phi$ can be defined as follows:
# 
# $\Phi^0_1$ corresponds to $\ \ \theta=0,\phi\in[0,2\pi],\ \ $ so
# $\ \ \Phi^0_1: \phi\to (0,0,1)$,
# 
# $\Phi^1_1$ corresponds to  $\ \ \theta=\frac{\pi}{2},\phi\in[0,2\pi], \ \ $ so $\Phi^1_1: \phi\to (r\cos\phi,r\sin\phi,0)$,
# 
# $\Phi^0_2$ corresponds to  $\ \ \theta\in[0,\frac{\pi}{2}],\phi=0, \ \ $ so $\Phi^0_2: \theta\to (r\sin\theta,0,r\cos\theta)$,
# 
# $\Phi^1_2$ corresponds to  $\ \ \theta\in[0,\frac{\pi}{2}],\phi=2/pi, \ \ $ so $\Phi^1_2: \theta\to (r\sin\theta,0,r\cos\theta)$.
# 
# Since by (20.16) $\int_{\partial\Phi}=\int_{\Phi^1_1}-\int_{\Phi^0_1} -\int_{\Phi^1_2}+\int_{\Phi^0_2}$, the  integral $\int_{\phi^0_1} a$ is zero and last two integrals cancel each other, the boundary integral reduces to $\int_{\Phi^1_1}a$. 
# 

# Define the pullback of $a$
# under $\psi:t\to (r\cos(t),\sin(t),0)$,

# In[20]:


N = Manifold(1, 'N')           # 1-dim manifold -curve
c.<t>=N.chart()                # polar coordinate on N
r=var('r')                               
                               # embedding N->R^3:
psi = N.diff_map(M, (r*cos(t), r*sin(t),0),
                 name='psi',latex_name=r'\psi')
psi.pullback(a).disp()         # pullback of a


# and use (reparametrized version of) (20.11), to obtain the boundary integral:
# $r^2\int_0^{2\pi}(\cos^2 t-\cos t\sin t)\;dt$.

# In[21]:


r^2*integrate(cos(t)^2-cos(t)*sin(t),(t,0,2*pi))


# Now we use the map 
# $\Phi:(\theta,\phi)\to 
# (r\sin(\theta)\cos(\phi), r\sin(\theta)\sin(\phi),r\cos(\theta))$,
# 

# In[22]:


r=var('r')
N = Manifold(2, 'N')              # manifold N=S^2
c_sph.<theta,phi>=N.chart()       # spherical coordinates
                                  # embedding S^2->R^3
Phi = N.diff_map(M, (r*sin(theta)*cos(phi), r*sin(theta)*sin(phi),r*cos(theta)),
                 name='psi',latex_name=r'\psi')


# and the pullback $\Phi^*da$.

# In[23]:


Phi.pullback(da).disp()


# We can apply (reparametrized version of) (20.11) to obtain the integral of $da$ over the hemisphere:
# $r^2\int_0^{2\pi}\int_0^{\pi/2}(\cos\phi\sin^2\phi+\cos\theta\sin\theta)\; d\theta d\phi.$

# In[24]:


r^2*integrate(integrate(cos(phi)*(sin(theta))^2+cos(theta)*sin(theta),
                        (theta,0,pi/2)),(phi,0,2*pi))


# <br>
# 
# **Example 20.9**
# 
# Compute the integral $\int_S x^2dy\wedge dz+\frac{1}{2}y^2dz\wedge dx+dx\wedge dy$ over the 
# tetrahedron bounded by the coordinate planes and the plane $x+y+z=1$.
# 
# <br>
# 
# As usual we start from defining 2-form $a$ and 3-form $da$.

# In[25]:


get_ipython().run_line_magic('display', 'latex')
M = Manifold(3, 'R^3')            # manifold M=R^3
c_xyz.<x,y,z> = M.chart()         # Cartesian coordinates
a=M.diff_form(2,name='a')
E=c_xyz.coframe()                 # coframe (dx,dy,dz)
dx,dy,dz=E[:]                     # define dx,dy,dz
a=x^2*dy.wedge(dz)+1/2*y^2*dz.wedge(dx)+dx.wedge(dy) #a
da=a.exterior_derivative()        # da
da.disp()


# Since the tetrahedron can be described by the inequalities
# 
# $$0\leq x\leq 1,\quad 0\leq y\leq 1-x,\quad 0\leq z\leq 1-x-y,$$
# 
# which is a  three-dimensional version of (20.25), by Stokes theorem, the boundary integral is equal to 
# 
# $$\int_0^1\Big[\int_0^{1-x}\Big[\int_0^{1-x-y} (2x+y)dz\Big] dy\Big] dx.$$

# In[26]:


integrate(integrate(integrate(2*x+y,(z,0,1-x-y)),(y,0,1-x)),(x,0,1))


# <br>
# 
# **Example 20.10**
# 
# Compute the integral $\int_S (x+y)dx+(2x-z)dy+(y+3z)dz,$ where $S$ is the boundary of the triangle with vertices $(3,0,0), \ (0,3,0), \ (0,0,3)$.
# 
# <br>
# 
# We start from defining 1-form $a$ and 2-form $da$.

# In[27]:


get_ipython().run_line_magic('display', 'latex')
M = Manifold(3, 'R^3')            # manifold M=R^3
c_xyz.<x,y,z> = M.chart()         # Cartesian coordinates
a=M.diff_form(2,name='a')
E=c_xyz.coframe()                 # coframe (dx,dy,dz)
dx,dy,dz=E[:]                     # define dx,dy,dz
a=(x+y)*dx+(2*x-z)*dy+(y+3*z)*dz  # a
da=a.exterior_derivative()        # da
da.disp()


# First we compute the integral of $da$ over the triangle using the pullback of $da$ under
# $\Phi:(X,Y)\to(X,Y,3-X-Y)$.

# In[28]:


N = Manifold(2, 'N')              # triangle S plane
c.<X,Y>=N.chart()                 # coord on N
                                  # N-> R^3
Phi = N.diff_map(M, (X, Y,3-X-Y),
                 name='Phi',latex_name=r'\Phi')
Phi.pullback(da).disp()           #pullback of da


# The iterated integral $\ \ \int_0^3[\int_0^{3-x}3dy]dx$ is equal to$\frac{27}{2}$.

# In[29]:


integrate(integrate(3,(y,0,3-x)),(x,0,3))


# The boundary of the triangle can be considered as the 1-chain $c=\psi_1+\psi_2+\psi_3$, where $\psi_1=-\Phi^1_2,\ \ \psi_2=-\Phi^0_1,\ \ \psi_3=\Phi^0_2$ are (reparametrized) singular 1-cubes (parametric curves).
# Note that the counterpart of $\Phi^1_1$ reduces to a point.
# 
# Let us define the space for the parameter $t$,

# In[30]:


N = Manifold(1, 'N')  # 1-dim manifold N "t axis" 
c.<t>=N.chart()       # coordinate on N


# next the singular 1-cubes $\psi_1,\psi_2,\psi_3$ and the pullbacks  $\psi_1^*a,\psi_2^*a,\psi_3^*a$:

# In[31]:


psi1 = N.diff_map(M, (3-t, t,0),
                 name='psi1',latex_name=r'\psi1')
psi1.pullback(a).disp()


# In[32]:


psi2 = N.diff_map(M, (0, 3-t,t),
                 name='psi2',latex_name=r'\psi2')
psi2.pullback(a).disp()


# In[33]:


psi3 = N.diff_map(M, (t, 0,3-t),
                 name='psi3',latex_name=r'\psi3')
psi3.pullback(a).disp()


# According to (reparametrized version of) (20.16), the integral of 1-form $a$ over 1-chain $c$ is equal to 
# 
# $$\textstyle\int_{\psi_1}a+\int_{\psi_2}a+\int_{\psi_3}a=\int_{[0,3]}\psi_1^*a+\int_{[0,3]}\psi_2^*a+\int_{[0,3]}\psi_3^*a\\
# =\int_0^3 ((-2t+3)+(3t+3)+(4t-9))dt.
# $$
# 

# In[34]:


integrate((-2*t+3)+(3*t+3)+(4*t-9),(t,0,3))


# <br>
# 
# ### Stokes theorem for chains
# 
# <br>
# 
# Extending Stokes theorem to $k$-chains is easy, since  for $k$-chain $c=\sum a_i\phi_i$, where $\phi_i$ are singular $k$-cubes in an open set $V$ and for  $(k-1)$-form $\omega$ in $V$,  we have by (20.14) and (20.23)
# 
# $$ \int_c d\omega=\sum a_i\int_{\phi_i} d\omega=\sum a_i\int_{\partial \phi_i}\omega=\int_{\partial c}\omega.$$

# <br>
# 
# **Remark.** In case of problems, some  iterated integrals can  be also computed with the help of  https://www.wolframalpha.com/ 
# 
# For example the command (without #):

# In[35]:


#integrate(integrate(integrate(2*exp(x+y)*sin(y), (z,0,1-x-y)), (y,0,1-x)), (x,0,1))


# in wolframalpha window returns
# 
# $$\int_0^1\int_0^{1-x}\int_0^{1-x-y}2\exp(x+y)\sin(y)\;dzdydx=
# -3-e(\sin(1)-2).$$

# **Remark**. Computing iterated integrals, SageMath (or strictly speaking Maxima) may need some assumptions, for example:

# In[36]:


get_ipython().run_line_magic('display', 'latex')
var('x y z')
forget()
assume(y+x-1<0)
assume(x-1<0)
integrate(integrate(integrate(2*exp(x+y+z)*sin(x-2*z),(z,0,1-x-y)),(y,0,1-x)),(x,0,1))


# ## What's next?
#   
# Take a look at the notebook [Connection](https://nbviewer.org/github/sagemanifolds/IntroToManifolds/blob/main/21Manifold_Connection.ipynb).