Skip to content

Instantly share code, notes, and snippets.

@icarob-eng
Created July 7, 2024 19:28
Show Gist options
  • Save icarob-eng/dc471edb1f04a5a43858da8c9153c597 to your computer and use it in GitHub Desktop.
Save icarob-eng/dc471edb1f04a5a43858da8c9153c597 to your computer and use it in GitHub Desktop.
integrais_multiplas.ipynb
Display the source blob
Display the rendered blob
Raw
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"provenance": [],
"authorship_tag": "ABX9TyP5YmRLsUwWOxglnTQQZL8E",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/icarob-eng/dc471edb1f04a5a43858da8c9153c597/integrais_multiplas.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "markdown",
"source": [
"# Cálculo númerico de expressões de integrais múltipllas em Python\n",
"\n",
"Este notebook iPy (`ipynb`) usa a biblioteca `sympy` para representar a integração simbolicamente e realiza a integração numericamente usando a biblioteca `scipy`. Para usar este notebook, você precisa saber pelo menos escrever expressões, usar funções, listas e, talvez classes, em Python. Ou não, você que sabe.\n",
"\n",
"Saiba que o `sympy` também tem seus próprios editores online com exemplos e vários tipos de saída interessantes em https://www.sympygamma.com/ e em https://live.sympy.org/.\n",
"\n",
"O trecho de código abaixo deve apenas ser executado, ele é responsável por importar bibliotecas e declarar a classe `IntegralIterada`.\n",
"\n",
"Crie uma cópia deste arquivo e vamos começar. Para executar uma célula, pode-se clicar no trângulo no canto superior esquerdo da tecla ou pressionar `ctrl + enter` tendo a célula selecionada.\n",
"\n",
"## Motivação do código\n",
"Este código seria inútil se fosse só para integrais simples, pois qualquer uma dessas coisas é fácil de se encontrar online, mas o motivo de eu tê-lo feito é porque é mais fácil de digitar expressões complexas de integrais iteradas, e se tem mais controle do ambiente de execução por aqui.\n",
"\n",
"Em geral eu fiz esse programa enquanto estudava, para saber se a mudança de variáveis ou ordem de integração estava certa, ou seja, produzia o mesmo resultado numérico, assim, o último exemplo é voltado a isto."
],
"metadata": {
"id": "9ERslQjdLJIM"
}
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"id": "Udyc0oJKwlHM"
},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.integrate import nquad\n",
"import sympy as sy\n",
"sy.init_printing()\n",
"\n",
"x, y, z = sy.symbols('x y z')\n",
"\n",
"class IntegralIterada:\n",
" def __init__(\n",
" self,\n",
" integrando: sy.Expr,\n",
" ordem: list[sy.Symbol],\n",
" lim_x: list[sy.Expr],\n",
" lim_y: list[sy.Expr]=None,\n",
" lim_z: list[sy.Expr]=None,\n",
" ):\n",
" \"\"\"\n",
" :param integrando: expressão a ser integrada\n",
" :param ordem: ordem de integração\n",
" :param lim_x: limites de x\n",
" :param lim_y: limites de y (opcional)\n",
" :param lim_z: limites de z (opcional)\n",
" \"\"\"\n",
" self.integrando = integrando\n",
" self.ordem = ordem\n",
" self.lims = {\n",
" x: lim_x,\n",
" y: lim_y,\n",
" z: lim_z,\n",
" }\n",
"\n",
" def expr(self):\n",
" return sy.Integral(\n",
" self.integrando,\n",
" *[(symbol, self.lims[symbol]) for symbol in self.ordem]\n",
" )\n",
"\n",
" def compute(self):\n",
" return nquad(\n",
" sy.lambdify(self.ordem, self.integrando),\n",
" [sy.lambdify(self.ordem[i+1:], self.lims[symbol])\n",
" for i, symbol in enumerate(self.ordem)]\n",
" )[0]"
]
},
{
"cell_type": "markdown",
"source": [
"## Criando expressões simbólicas\n",
"Para vizualizar as expressões matemáticas usadas, as escrevemos usando [`sympy`](https://www.sympy.org/), assim, se quisermos *representar* a raíz de um número, ou o seno de um número, usamos:"
],
"metadata": {
"id": "Al5AIEnFQ917"
}
},
{
"cell_type": "code",
"source": [
"sy.sqrt(2)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 38
},
"id": "y9MvokecPcJA",
"outputId": "bd6941be-053e-422e-cdc0-b9385b240fc6"
},
"execution_count": 2,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"√2"
],
"text/latex": "$\\displaystyle \\sqrt{2}$"
},
"metadata": {},
"execution_count": 2
}
]
},
{
"cell_type": "code",
"source": [
"sy.sin(2)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 37
},
"id": "k4Zx8mSSSMoq",
"outputId": "0122baf1-2671-4c77-8592-aa4117b99086"
},
"execution_count": 3,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"sin(2)"
],
"text/latex": "$\\displaystyle \\sin{\\left(2 \\right)}$"
},
"metadata": {},
"execution_count": 3
}
]
},
{
"cell_type": "markdown",
"source": [
"O `sympy` possui a representação simbólica de diversas funções e constantes, como mostraremos a seguir. Quais? Muitas, consulte a documentação (https://docs.sympy.org/latest/index.html). Veja que o `sympy`, por padrão, executa a computação de alguns valores inteiros, mas deixa racionais em forma simbólica ao invés de computar o número."
],
"metadata": {
"id": "-6U8aUgXSUlw"
}
},
{
"cell_type": "code",
"source": [
"sy.root(2, 3) # resultado simbólico"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 38
},
"id": "TfRcojvsS5_W",
"outputId": "a4732752-e269-4d2e-a99f-1757374da3a2"
},
"execution_count": 4,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"3 ___\n",
"╲╱ 2 "
],
"text/latex": "$\\displaystyle \\sqrt[3]{2}$"
},
"metadata": {},
"execution_count": 4
}
]
},
{
"cell_type": "code",
"source": [
"sy.root(8, 3) # resultado numérico"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 37
},
"id": "XRq2eZZMS8_A",
"outputId": "a22d4ff0-ddea-471c-ee03-a18f48a53e4c"
},
"execution_count": 5,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2"
],
"text/latex": "$\\displaystyle 2$"
},
"metadata": {},
"execution_count": 5
}
]
},
{
"cell_type": "code",
"source": [
"sy.pi ** 2 # resultado simbólico"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 38
},
"id": "20-ponHxTFte",
"outputId": "7814d6cf-9150-400a-8003-f23c0b0152f0"
},
"execution_count": 6,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
" 2\n",
"π "
],
"text/latex": "$\\displaystyle \\pi^{2}$"
},
"metadata": {},
"execution_count": 6
}
]
},
{
"cell_type": "code",
"source": [
"sy.tan(sy.pi) # resultado numérico"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 37
},
"id": "R2UZMRjQTOyO",
"outputId": "669d7c3f-a245-473c-add9-e0b4c943ce0f"
},
"execution_count": 7,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"0"
],
"text/latex": "$\\displaystyle 0$"
},
"metadata": {},
"execution_count": 7
}
]
},
{
"cell_type": "markdown",
"source": [
"Se você quiser saber o valor numérico do resultado, basta acrescentar `.doit()` (faça-o) ao final da expressão."
],
"metadata": {
"id": "PEJPEjPRTbdF"
}
},
{
"cell_type": "code",
"source": [
"(sy.pi ** 2).doit()"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 38
},
"id": "NqznFzIJTqMA",
"outputId": "33f713eb-ab70-473f-d36c-417f15e152dd"
},
"execution_count": 8,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
" 2\n",
"π "
],
"text/latex": "$\\displaystyle \\pi^{2}$"
},
"metadata": {},
"execution_count": 8
}
]
},
{
"cell_type": "markdown",
"source": [
"Como tudo no Python, é possível atribuir uma variável *do Python* a uma expressão:"
],
"metadata": {
"id": "yAsALnBHVqh7"
}
},
{
"cell_type": "code",
"source": [
"a = sy.Rational(2, 3)"
],
"metadata": {
"id": "-dQWtnEMVxpI"
},
"execution_count": 9,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"Veja que quando se atribui uma expressão a uma variável, o notebook não mostra a expressão. Para mostrar uma expressão, basta deixar a variável (ou mesmo a expressão) na última linha de uma célula."
],
"metadata": {
"id": "vEXvI6iFWC7f"
}
},
{
"cell_type": "code",
"source": [
"a = sy.Rational(2, 3)\n",
"a"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 51
},
"id": "sDQvWMRcW9io",
"outputId": "b1c797a5-e394-423c-9ff2-e3a547816c45"
},
"execution_count": 10,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"2/3"
],
"text/latex": "$\\displaystyle \\frac{2}{3}$"
},
"metadata": {},
"execution_count": 10
}
]
},
{
"cell_type": "markdown",
"source": [
"Note que para reprentar uma fração entre dois números (não simbólicos), você deve usar `sy.Rational(2, 3)`. Cuidado para não se confundir e digitar `sy.Rational(2/3)`, pois, por limitações computacionais, isto gera um resultado diferente, visto que o primeiro é uma representação simbólica da fração e o segundo um valor real com precisão limitada no computador. E apenas `2/3` também. Sinta-se a vontade para testar."
],
"metadata": {
"id": "Lzqu74u_XPBi"
}
},
{
"cell_type": "markdown",
"source": [
"### Variáveis simbólicas\n",
"No presente notebook, já criamos as variáveis *simbólicas* $x$, $y$ e $z$ e atribuímos às variáveis do Python `x`, `y` e `z`, por conveniência e também porque que são usadas como possíveis variáveis de integração na classe `IntegralIterada`. Então, se pretende usar coorddenadas polares, por favor escreva elas em termos de xyz, por exemplo:\n",
"\n",
"|rec|pol|\n",
"|---|---|\n",
"|$x$|$r$|\n",
"|$y$|$θ$|\n",
"|$z$|$ϕ$|\n",
"\n",
"> Nota: *variável simbólica* é uma classe do `sympy` (especificamente `Symbol`) que por padrão não possue valor numérico e só pode ser suada em expressões. Não confuda com *variável do Python* que representa um determinado espaço na memória do computador que armazena qualquer tipo de dado **incluindo variaveis simbólicas**.\n",
"\n",
"Com as variáveis é possível escrever expressões sem valor definido (um lógico chamaria de \"predicado\"):"
],
"metadata": {
"id": "8kQwYzz6Txxa"
}
},
{
"cell_type": "code",
"source": [
"2*(3*x+1)*sy.sqrt(4*y)*3"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 38
},
"id": "_U-R86gkU-9A",
"outputId": "558bc792-d1a2-4b34-85da-bed99856a891"
},
"execution_count": 11,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6⋅√y⋅(6⋅x + 2)"
],
"text/latex": "$\\displaystyle 6 \\sqrt{y} \\left(6 x + 2\\right)$"
},
"metadata": {},
"execution_count": 11
}
]
},
{
"cell_type": "code",
"source": [
"(2*(3*x+1)*sy.sqrt(4*y)*3).doit() # não tem valor numérico"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 38
},
"id": "M5TQC5nSVDhH",
"outputId": "8898c390-c4e7-4664-c8f6-5914cfa726e6"
},
"execution_count": 12,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"6⋅√y⋅(6⋅x + 2)"
],
"text/latex": "$\\displaystyle 6 \\sqrt{y} \\left(6 x + 2\\right)$"
},
"metadata": {},
"execution_count": 12
}
]
},
{
"cell_type": "markdown",
"source": [
"## Integração simbólica e numérica\n",
"O `sympy` permite representar integrais com facilidade, mas computa-las, seja numéricamente, e ainda mais simbolicamente, é um pouco mais difícil. Por isto criei a classe `IntegralIterada`, que converte internamente a expressão simbólica em uma função do Python e a integra numericamente. Como tudo no Python, você pode atribuir uma instância de `IntegralIterada` a uma variável e fazer operações com a integral em si."
],
"metadata": {
"id": "-053cA67ZHm9"
}
},
{
"cell_type": "markdown",
"source": [
"### Integral simples\n",
"Comecemos com um exemplo: na célula abaixo criamos uma instância de integral e atribuimos à variável `a`; em seguida, apresentamos a expressão que foi passada, com os limites de integração e variáveis de integração aplicada. Isto é útil para confirmar o que estamos escrevendo."
],
"metadata": {
"id": "TQpGtfT4aPiO"
}
},
{
"cell_type": "code",
"source": [
"a = IntegralIterada(\n",
" integrando=4*x*sy.sqrt(3*x+1)*4 + 2*(3*x+1)*sy.sqrt(4*x)*3,\n",
" ordem=[x],\n",
" lim_x=[0, 1]\n",
")\n",
"a.expr()"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 80
},
"id": "QjccR-v-Me0W",
"outputId": "449ad4f3-84ff-4621-fc5e-bee1d7dcbe19"
},
"execution_count": 13,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"1 \n",
"⌠ \n",
"⎮ ⎛ _________⎞ \n",
"⎮ ⎝6⋅√x⋅(6⋅x + 2) + 16⋅x⋅╲╱ 3⋅x + 1 ⎠ dx\n",
"⌡ \n",
"0 "
],
"text/latex": "$\\displaystyle \\int\\limits_{0}^{1} \\left(6 \\sqrt{x} \\left(6 x + 2\\right) + 16 x \\sqrt{3 x + 1}\\right)\\, dx$"
},
"metadata": {},
"execution_count": 13
}
]
},
{
"cell_type": "markdown",
"source": [
"Veja que os parâmetros da função `IntegralIterada` são:\n",
"- `integrando`, que deve ser uma expressão simbólica;\n",
"- `ordem` de integração: a lista ordenada de variáveis de integração, que devem ser as variáveis `x`, `y` ou `z`;\n",
"- E `lim_x` que é um par com os limites de integração;\n",
"\n",
"O outro recurso dessa classe é `.compute()`, para obter o valor numérico da integral. Isto é utíl para verificar se duas integrais estão numericamente equivalentes considerando os limites, mas saiba que o `sympy` permite comparar expressões, (com algumas limitações"
],
"metadata": {
"id": "0dnD6KdQbZlm"
}
},
{
"cell_type": "code",
"source": [
"a.compute()"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 37
},
"id": "K234TYzJSkFQ",
"outputId": "736156ca-51f3-4661-976a-bebb529effee"
},
"execution_count": 14,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"36.14814814814814"
],
"text/latex": "$\\displaystyle 36.1481481481481$"
},
"metadata": {},
"execution_count": 14
}
]
},
{
"cell_type": "markdown",
"source": [
"### Integrais múltiplas\n",
"Vamos testar uma mudança de ordem de integração da seguinte integral dupla. Como é uma mudança de ordem de integração, podemos armazenar a expressão do integrando numa variável para usar nas duas integrais e diminuira chance de erro de digitação."
],
"metadata": {
"id": "MBasgaSefCl-"
}
},
{
"cell_type": "code",
"source": [
"expr = sy.cos(x) * sy.sqrt(1 + sy.cos(x)**2)\n",
"expr"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 40
},
"id": "0p190CmRi-B7",
"outputId": "40730fa2-f598-4c78-e401-f3b8a7403491"
},
"execution_count": 15,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
" _____________ \n",
" ╱ 2 \n",
"╲╱ cos (x) + 1 ⋅cos(x)"
],
"text/latex": "$\\displaystyle \\sqrt{\\cos^{2}{\\left(x \\right)} + 1} \\cos{\\left(x \\right)}$"
},
"metadata": {},
"execution_count": 15
}
]
},
{
"cell_type": "code",
"source": [
"a = IntegralIterada(\n",
" integrando=expr,\n",
" ordem=[x, y],\n",
" lim_x=[sy.asin(y), sy.pi/2],\n",
" lim_y=[0, 1]\n",
")\n",
"a.expr()"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 87
},
"id": "NeYSxFJYNq-f",
"outputId": "b1814830-bf68-498c-9658-33a2d47a6435"
},
"execution_count": 16,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
" π \n",
" ─ \n",
"1 2 \n",
"⌠ ⌠ \n",
"⎮ ⎮ _____________ \n",
"⎮ ⎮ ╱ 2 \n",
"⎮ ⎮ ╲╱ cos (x) + 1 ⋅cos(x) dx dy\n",
"⌡ ⌡ \n",
"0 asin(y) "
],
"text/latex": "$\\displaystyle \\int\\limits_{0}^{1}\\int\\limits_{\\operatorname{asin}{\\left(y \\right)}}^{\\frac{\\pi}{2}} \\sqrt{\\cos^{2}{\\left(x \\right)} + 1} \\cos{\\left(x \\right)}\\, dx\\, dy$"
},
"metadata": {},
"execution_count": 16
}
]
},
{
"cell_type": "markdown",
"source": [
"Agora, queremos testar se a seguinte integral dupla é equivalente à mudança de variável da integral acima."
],
"metadata": {
"id": "H_MXZ9Efijv4"
}
},
{
"cell_type": "code",
"source": [
"b = IntegralIterada(\n",
" integrando=expr,\n",
" ordem=[y, x],\n",
" lim_x=[0, sy.pi/2],\n",
" lim_y=[0, sy.sin(x)]\n",
")\n",
"b.expr()"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 83
},
"id": "U6BEjE1PiRmL",
"outputId": "2809de77-bc14-4f80-9d47-e8a23b90e004"
},
"execution_count": 17,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"π \n",
"─ \n",
"2 sin(x) \n",
"⌠ ⌠ \n",
"⎮ ⎮ _____________ \n",
"⎮ ⎮ ╱ 2 \n",
"⎮ ⎮ ╲╱ cos (x) + 1 ⋅cos(x) dy dx\n",
"⌡ ⌡ \n",
"0 0 "
],
"text/latex": "$\\displaystyle \\int\\limits_{0}^{\\frac{\\pi}{2}}\\int\\limits_{0}^{\\sin{\\left(x \\right)}} \\sqrt{\\cos^{2}{\\left(x \\right)} + 1} \\cos{\\left(x \\right)}\\, dy\\, dx$"
},
"metadata": {},
"execution_count": 17
}
]
},
{
"cell_type": "markdown",
"source": [
"Para saber se os limites são equivalentes, podemos avaliar as expressões numericamente e comparar a diferença entre os dois (o resíduo):"
],
"metadata": {
"id": "vuYjFe-DjLbd"
}
},
{
"cell_type": "code",
"source": [
"print(f'''Integral a: {a.compute()}\n",
" \\nIntegral b: {b.compute()}\n",
" \\nResíduo: {a.compute() - b.compute():5f}''')"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "5xjgzD0jWXRI",
"outputId": "7d57371f-b9d0-4bf6-e08b-3adeb41f10d1"
},
"execution_count": 18,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Integral a: 0.60947570824873\n",
" \n",
"Integral b: 0.60947570824873\n",
" \n",
"Resíduo: 0.000000\n"
]
}
]
},
{
"cell_type": "markdown",
"source": [
"# Sua vez\n",
"\n",
"Pode fazer seus estudos em qualquer parte do notebook, ou então usar o espaço abaixo para fazer os próprios testes."
],
"metadata": {
"id": "C3ouZIFwjiq1"
}
},
{
"cell_type": "code",
"source": [],
"metadata": {
"id": "hEQmHWSVjmt5"
},
"execution_count": 18,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment