# SymPyLab

**INTEGRANTES**

Ivonne Nathalia Uribe Lara

David Antonio Garzón Avendaño

Juan Carlos Sanchez Orjuela


SymPy’s documentation
- https://docs.sympy.org/latest/index.html


## SymPy’s polynomials 
- https://docs.sympy.org/latest/modules/polys/basics.html#polynomials 

- (x-1)(x-2)(x-3)(x-4)(x-5) = x^5 - 15 x^4  + 85 x^3 - 225 x^2 + 274 x - 120

- (x^5 - 15 x^4  + 85 x^3 - 225 x^2 + 274 x - 120) / (x-1) = x^4  - 14 x^3  + 71 x^2  - 154 x + 120

 



```
from sympy import Symbol
from sympy import div

x = Symbol('x')

p = x**5 - 15*x**4  + 85*x**3 - 225*x**2 + 274*x - 120

p, r = div(p,  x-1)

print(p)
print(r)

p, r = div(p,  x-2)

print(p)
print(r)

p, r = div(p,  x-3)

print(p)
print(r)

p, r = div(p,  x-4)

print(p)
print(r)


```

    x**4 - 14*x**3 + 71*x**2 - 154*x + 120
    0
    x**3 - 12*x**2 + 47*x - 60
    0
    x**2 - 9*x + 20
    0
    x - 5
    0


Edición

 Polinomio:
 
 (x-1)(x-2)(x-3)(x-4)(x-5)(x-6)(x-7)(x-8)(x-9)(x-10) = x^10 - 55x^9 = 1320x^8 - 18150x^7 + 157773x^6 - 902055x^5 + 3416930x^4 + 8409500x^3 - 10628640x + 3628800


(x-1)(x-2)(x-3)(x-4)(x-5)(x-6)(x-7)(x-8)(x-9)(x-10)/(x-1) = x^9 - 54x^8 + 1266x^7 - 16884x^6 + 140889x^5 - 761166x^4 + 2655764^x^3 - 5753736^x^2 + 6999840^x - 3628800


wolfram:
 


 


```
#edicion

from sympy import Symbol
from sympy import div

x = Symbol('x')

p = (x-1)*(x-2)*(x-3)*(x-4)*(x-5)*(x-6)*(x-7)*(x-8)*(x-9)*(x-10)

p, r = div(p,  x-1)

print(p)
print(r)

p, r = div(p,  x-2)

print(p)
print(r)

p, r = div(p,  x-3)

print(p)
print(r)

p, r = div(p,  x-4)

print(p)
print(r)

p, r = div(p,  x-5)

print(p)
print(r)

p, r = div(p,  x-6)

print(p)
print(r)

p, r = div(p,  x-7)

print(p)
print(r)

p, r = div(p,  x-8)

print(p)
print(r)

p, r = div(p,  x-9)

print(p)
print(r)


```

    x**9 - 54*x**8 + 1266*x**7 - 16884*x**6 + 140889*x**5 - 761166*x**4 + 2655764*x**3 - 5753736*x**2 + 6999840*x - 3628800
    0
    x**8 - 52*x**7 + 1162*x**6 - 14560*x**5 + 111769*x**4 - 537628*x**3 + 1580508*x**2 - 2592720*x + 1814400
    0
    x**7 - 49*x**6 + 1015*x**5 - 11515*x**4 + 77224*x**3 - 305956*x**2 + 662640*x - 604800
    0
    x**6 - 45*x**5 + 835*x**4 - 8175*x**3 + 44524*x**2 - 127860*x + 151200
    0
    x**5 - 40*x**4 + 635*x**3 - 5000*x**2 + 19524*x - 30240
    0
    x**4 - 34*x**3 + 431*x**2 - 2414*x + 5040
    0
    x**3 - 27*x**2 + 242*x - 720
    0
    x**2 - 19*x + 90
    0
    x - 10
    0


## SymPy’s polynomial simple univariate polynomial factorization
- https://docs.sympy.org/latest/modules/polys/wester.html#simple-univariate-polynomial-factorization
- factor(x\*\*5 - 15\*x\*\*4  + 85\*x\*\*3 - 225\*x\*\*2 + 274\*x - 120)





```
from sympy import *
x = Symbol('x')
factor(x**5 - 15*x**4  + 85*x**3 - 225*x**2 + 274*x - 120)
```




    (x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1)



Edición

factor(x^4 - 46x^3 + 233x^2 + 5656x + 9996)





```
from sympy import *
x = Symbol('x')
factor(x**4 - 46*x**3  +233*x**2 + 5656*x + 9996)
```




    (x - 34)*(x - 21)*(x + 2)*(x + 7)



## SymPy’s solvers
- https://docs.sympy.org/latest/tutorial/solvers.html
- x\*\*5 - 15\*x\*\*4  + 85\*x\*\* 3 - 225\*x\*\* 2 + 274\*x - 120 = 0




```
from sympy import *
x = Symbol('x')
solveset(Eq(x**5 - 15*x**4  + 85*x**3 - 225*x**2 + 274*x - 120, 0), x)

```




    {1, 2, 3, 4, 5}



Edición 


resolve(x^4 - e^x - pi = 0)

Wolfram:




```
#edicion 

from sympy import *
import numpy as np
x = Symbol('x')
solveset(Eq(x**4 -32*x**3 - np.pi , 0), x)
```




    {-0.459132099076072, 32.0000958729375, 0.229518113069273 - 0.401431993025616*I, 0.229518113069273 + 0.401431993025616*I}



## SymPy’s Symbolic and Numercical Complex Evaluations
- https://docs.sympy.org/latest/modules/evalf.html](https://)
- x = x1 + I*x2,y = y1 + I*y2, z = z1 + I*z2, x*y*z




```
from sympy import *
x1, x2, y1, y2, z1, z2 = symbols("x1 x2 y1 y2 z1 z2", real=True)  
x = x1 + I*x2
y = y1 + I*y2
z = z1 + I*z2

print(x*y*z)
print(expand(x*y*z))
print(expand((x*y)*z))
print(expand(x*(y*z)))

w=N(1/(pi + I), 20)
print('w=',w)
```

    (x1 + I*x2)*(y1 + I*y2)*(z1 + I*z2)
    x1*y1*z1 + I*x1*y1*z2 + I*x1*y2*z1 - x1*y2*z2 + I*x2*y1*z1 - x2*y1*z2 - x2*y2*z1 - I*x2*y2*z2
    x1*y1*z1 + I*x1*y1*z2 + I*x1*y2*z1 - x1*y2*z2 + I*x2*y1*z1 - x2*y1*z2 - x2*y2*z1 - I*x2*y2*z2
    x1*y1*z1 + I*x1*y1*z2 + I*x1*y2*z1 - x1*y2*z2 + I*x2*y1*z1 - x2*y1*z2 - x2*y2*z1 - I*x2*y2*z2
    w= 0.28902548222223624241 - 0.091999668350375232456*I


Edición 

x = x1 + I*x2, y = y1 + I*y2, z = z1 + I*z2, a = a1 + I*a2, x*y*z/a

Wolfram





```
from sympy import *
x1, x2, y1, y2, z1, z2, a1, a2 = symbols("x1 x2 y1 y2 z1 z2 a1 a2", real=True)  
x = x1 + I*x2
y = y1 + I*y2
z = z1 + I*z2
a = a1 + I*a2

print(x*y*(z/a))
print(expand(x*y*(z/a)))
print(expand((x*y)*(z/a)))
print(expand((z/a)*x*y))


wres =N(pi**2 - sin(2), 128)
print('Wres=',wres)

```

    (x1 + I*x2)*(y1 + I*y2)*(z1 + I*z2)/(a1 + I*a2)
    x1*y1*z1/(a1 + I*a2) + I*x1*y1*z2/(a1 + I*a2) + I*x1*y2*z1/(a1 + I*a2) - x1*y2*z2/(a1 + I*a2) + I*x2*y1*z1/(a1 + I*a2) - x2*y1*z2/(a1 + I*a2) - x2*y2*z1/(a1 + I*a2) - I*x2*y2*z2/(a1 + I*a2)
    x1*y1*z1/(a1 + I*a2) + I*x1*y1*z2/(a1 + I*a2) + I*x1*y2*z1/(a1 + I*a2) - x1*y2*z2/(a1 + I*a2) + I*x2*y1*z1/(a1 + I*a2) - x2*y1*z2/(a1 + I*a2) - x2*y2*z1/(a1 + I*a2) - I*x2*y2*z2/(a1 + I*a2)
    x1*y1*z1/(a1 + I*a2) + I*x1*y1*z2/(a1 + I*a2) + I*x1*y2*z1/(a1 + I*a2) - x1*y2*z2/(a1 + I*a2) + I*x2*y1*z1/(a1 + I*a2) - x2*y1*z2/(a1 + I*a2) - x2*y2*z1/(a1 + I*a2) - I*x2*y2*z2/(a1 + I*a2)
    Wres= 8.9603069742636769234384711339644062926114444357929003580343763646890775208784216983816467147935928515140580170316236578164913950


## SymPy’s integrals
- https://docs.sympy.org/latest/modules/integrals/integrals.html
- [risk-engineering.org](https://risk-engineering.org/notebook/monte-carlo-LHS.html)



Let’s start with a simple integration problem in 1D,

$$\int_1^5 x^2 dx$$
 
This is easy to solve analytically, and we can use the SymPy library in case you’ve forgotten how to resolve simple integrals.


```
import sympy
# we’ll save results using different methods in this data structure, called a dictionary
result = {}  
x = sympy.Symbol("x")
i = sympy.integrate(x**2)
print(i)
result["analytical"] = float(i.subs(x, 5) - i.subs(x, 1))
print("Analytical result: {}".format(result["analytical"]))
```

    x**3/3
    Analytical result: 41.333333333333336


Edición

Inetgral: 

$$\int_3^8 x^3 + 8 dx$$


```
#Edicion

import sympy
resultado = {}  
x = sympy.Symbol("x")
i = sympy.integrate(x**3 + 8)
print(i)
resultado["ant"] = float(i.subs(x, 8) - i.subs(x, 3))
print("Resiltado: {}".format(resultado["ant"]))
```

    x**4/4 + 8*x
    Resiltado: 1043.75


**Integrating with Monte Carlo** 
[risk-engineering.org](https://risk-engineering.org/notebook/monte-carlo-LHS.html) 

We can estimate this integral using a standard Monte Carlo method, where we use the fact that the expectation of a random variable is related to its integral

$$\mathbb{E}(f(x)) = \int_I f(x) dx $$

We will sample a large number N of points in I and calculate their average, and multiply by the range over which we are integrating.


```
import numpy
N = 1_000_000
accum = 0
for i in range(N):
    x = numpy.random.uniform(1, 5)
    accum += x**2
volume = 5 - 1
result["MC"] = volume * accum / float(N)
print("Standard Monte Carlo result: {}".format(result["MC"]))
```

    Standard Monte Carlo result: 41.32998559988266


Edición

Estimamos la siguiente integral 

$$\int_3^8 x^3 + 8 dx$$


por Monte Carlo 


```
import numpy 
N = 10_000_000
temp = 0
for i in range(N):
    x = numpy.random.uniform(3, 8)
    temp += x**3 + 8
volume = 8 - 3
resultado["MontCar"] = volume * temp / float(N)
print("Resultado por Monte Carlo: {}".format(resultado["MontCar"]))
```

    Resultado por Monte Carlo: 1043.9422026850486


- integrate(x\*\*2 * sin(x)\*\*3)




```
import sympy
x = Symbol("x")
i = integrate(x**2 * sin(x)**3)
print(i)
print(float(i.subs(x, 5) - i.subs(x, 1)))
```

    -x**2*sin(x)**2*cos(x) - 2*x**2*cos(x)**3/3 + 14*x*sin(x)**3/9 + 4*x*sin(x)*cos(x)**2/3 + 14*sin(x)**2*cos(x)/9 + 40*cos(x)**3/27
    -15.42978215330555


Edición

Inetgramos : sin(x) + (cos(x)^2)


Wolfram: 




```
#Edicion 
x = Symbol("x")
i = integrate(sin(x) + (cos(x)**2))
print(i)
print(float(i.subs(x, 8) - i.subs(x, 3)))

```

    x/2 + sin(x)*cos(x)/2 - cos(x)
    1.6533855825916333



```
import numpy
N = 100_000
accum = 0
l =[]
for i in range(N):
    x = numpy.random.uniform(1, 5)
    accum += x**2 * sin(x)**3
volume = 5 - 1
result["MC"] = volume * accum / float(N)
print("Standard Monte Carlo result: {}".format(result["MC"]))
```

    Standard Monte Carlo result: -15.3350781160391



```
#Edicion 

import numpy 
N = 1_000_000
temp = 0
for i in range(N):
    x = numpy.random.uniform(3, 8)
    temp += numpy.sin(x) + (numpy.cos(x))**2
volume = 8 - 3
resultado["MontCar2"] = volume * temp / float(N)
print("Resultado por Monte Carlo: {}".format(resultado["MontCar2"]))
```

    Resultado por Monte Carlo: 1.644611681543184


**A higher dimensional integral** [risk-engineering.org](https://risk-engineering.org/notebook/monte-carlo-LHS.html) 


Let us now analyze an integration problem in dimension 4, the Ishigami function. This is a well-known function in numerical optimization and stochastic analysis, because it is very highly non-linear.


```
import sympy

x1 = sympy.Symbol("x1")
x2 = sympy.Symbol("x2")
x3 = sympy.Symbol("x3")
expr = sympy.sin(x1) + 7*sympy.sin(x2)**2 + 0.1 * x3**4 * sympy.sin(x1)
res = sympy.integrate(expr,
                      (x1, -sympy.pi, sympy.pi),
                      (x2, -sympy.pi, sympy.pi),
                      (x3, -sympy.pi, sympy.pi))
# Note: we use float(res) to convert res from symbolic form to floating point form
result = {} 
result["analytical"] = float(res)
print("Analytical result: {}".format(result["analytical"]))
```

    Analytical result: 868.175747048395



```
#Edicion 

import sympy 

x = sympy.Symbol("x") 
y = sympy.Symbol("y") 
z = sympy.Symbol("z")
expr = sympy.sin(x) + 8*sympy.sin(y)**2 + 0.34 * z**4 * sympy.sin(x) 
res = sympy.integrate(expr,
                      (x, -sympy.pi, sympy.pi),
                      (y, -sympy.pi, sympy.pi),
                      (z, -sympy.pi, sympy.pi))

resultados = {} 
resultados["analytical"] = float(res)
print("Analytical result: {}".format(resultados["analytical"]))

```

    Analytical result: 992.2008537695942



```
N = 10_000
accum = 0
for i in range(N):
    xx1 = numpy.random.uniform(-numpy.pi, numpy.pi)
    xx2 = numpy.random.uniform(-numpy.pi, numpy.pi)
    xx3 = numpy.random.uniform(-numpy.pi, numpy.pi)
    accum += numpy.sin(xx1) + 7*numpy.sin(xx2)**2 + 0.1 * xx3**4 * numpy.sin(xx1)
volume = (2 * numpy.pi)**3
result = {} 
result["MC"] = volume * accum / float(N)
print("Standard Monte Carlo result: {}".format(result["MC"]))
```

    Standard Monte Carlo result: 860.6796782529639



```
#Edicion 

import numpy 
N = 100_000
temp = 0
for i in range(N):
  xx = numpy.random.uniform(-numpy.pi, numpy.pi)
  yy = numpy.random.uniform(-numpy.pi, numpy.pi)
  zz = numpy.random.uniform(-numpy.pi, numpy.pi)
  temp += np.sin(xx) + 8*np.sin(yy)**2 + 0.34 * zz**4 * np.sin(xx)

volumen = (2 * numpy.pi)**3
resultados["MonteCarlo"] = volumen * temp / float(N)
print("Resultado Usando Montecarlo: {}".format(resultados["MonteCarlo"]))

```

    Resultado Usando Montecarlo: 983.0796811421585



```
import math
import numpy
# adapted from https://mail.scipy.org/pipermail/scipy-user/2013-June/034744.html
def halton(dim: int, nbpts: int):
    h = numpy.full(nbpts * dim, numpy.nan)
    p = numpy.full(nbpts, numpy.nan)
    P = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
    lognbpts = math.log(nbpts + 1)
    for i in range(dim):
        b = P[i]
        n = int(math.ceil(lognbpts / math.log(b)))
        for t in range(n):
            p[t] = pow(b, -(t + 1))

        for j in range(nbpts):
            d = j + 1
            sum_ = math.fmod(d, b) * p[0]
            for t in range(1, n):
                d = math.floor(d / b)
                sum_ += math.fmod(d, b) * p[t]

            h[j*dim + i] = sum_
    return h.reshape(nbpts, dim)
```


```
import matplotlib.pyplot as plt
N = 1000
seq = halton(2, N)
plt.title("2D Halton sequence")
# Note: we use "alpha=0.5" in the scatterplot so that the plotted points are semi-transparent
# (alpha-transparency of 0.5 out of 1), so that we can see when any points are superimposed.
plt.axes().set_aspect('equal')
plt.scatter(seq[:,0], seq[:,1], marker=".", alpha=0.5);
```


```
#Edicion
import numpy
import math 

def halton2(dimencion: int, nbpts: int):
    h = numpy.full(nbpts * dimencion, numpy.nan)
    p = numpy.full(nbpts, numpy.nan)
    P = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    lognbpts = math.log(nbpts + 1)
    for i in range(dimencion):
        b = P[i]
        n = int(math.ceil(lognbpts / math.log(b)))
        for t in range(n):
            p[t] = b**(-(t + 1))

        for j in range(nbpts):
            d = j + 1
            sum_ = math.fmod(d, b) * p[0]
            for t in range(1, n):
                d = math.floor(d / b)
                sum_ += math.fmod(d, b) * p[t]

            h[j*dimencion + i] = sum_
    return h.reshape(nbpts, dimencion)

```


```
#Edicion 

import matplotlib.pyplot as plt
N = 1000
seq = halton2(2, N)
plt.title("Secuencia de Halton")
plt.axes().set_aspect('equal')
plt.scatter(seq[:,0], seq[:,1], marker=".", alpha=0.5, color="green");
```


```
N = 10_000

seq = halton(3, N)
accum = 0
for i in range(N):
    xx1 = -numpy.pi + seq[i][0] * numpy.pi * 2
    xx2 = -numpy.pi + seq[i][1] * numpy.pi * 2
    xx3 = -numpy.pi + seq[i][2] * numpy.pi * 2
    accum += numpy.sin(xx1) + 7*numpy.sin(xx2)**2 + 0.1 * xx3**4 * numpy.sin(xx1)
volume = (2 * numpy.pi)**3
result = {} 
result["MC"] = volume * accum / float(N)
print("Qausi Monte Carlo Halton Sequence result: {}".format(result["MC"]))
```

    Qausi Monte Carlo Halton Sequence result: 868.238928030592



```
#Edicion 

import numpy as np
N = 100_000
temp = 0
seq2 = halton2(3,N)
for i in range(N):
  xx = -numpy.pi + seq2[i][0] * numpy.pi * 2
  yy = -numpy.pi + seq2[i][1] * numpy.pi * 2
  zz = -numpy.pi + seq2[i][2] * numpy.pi * 2
  temp += np.sin(xx) + 8*np.sin(yy)**2 + 0.34 * zz**4 * np.sin(xx)

resultados ={}
volumen = (2 * numpy.pi)**3
resultados["MonteCarlo"] = volumen * temp / float(N)
print("Resultado Usando Montecarlo: {}".format(resultados["MonteCarlo"]))
```

    Resultado Usando Montecarlo: 992.1898731825128


Como punto adicional, vamos a hacer algunas integrales dobles y triples (no se anexan fotos de Wolfram debido a tiene problemas para resolver este tipo de integrales)


- $$\int_2^5 \int_3^6 x^2 + sin^2(y) dx dy$$

- $$\int_3^5 \int_6^7 ln(x) + cos(y) dx dy$$

- $$\int_1^5 \int_1^2 \int_2^3 e^x + ln(y) - cos(z) dx dy dz $$



```
#Edicion 

import sympy

x = sympy.Symbol("x")
y = sympy.Symbol("y")
expr = x**2 + sin(y)**2
res = sympy.integrate(expr,
                      (x, 2, 5),
                      (y, 3, 6))
result = {} 
result["analytical"] = float(res)
print("Analytical result: {}".format(result["analytical"]))
```

    Analytical result: 121.69286806485113



```
#Edicion 

import sympy

x = sympy.Symbol("x")
y = sympy.Symbol("y")
expr = ln(x) + cos(x)
res = sympy.integrate(expr,
                      (x, 3, 5),
                      (y, 6, 7))
result = {} 
result["analytical"] = float(res)
print("Analytical result: {}".format(result["analytical"]))
```

    Analytical result: 1.651308413443167



```
#Edicion 
import math
x = sympy.Symbol("x")
y = sympy.Symbol("y")
z = sympy.Symbol("z")
expr = math.e**x + ln(y) - cos(z)
res = sympy.integrate(expr,
                      (x, 1, 5),
                      (y, 1, 2),
                      (z, 2, 3))

result = {} 
result["analytical"] = float(res)
print("Analytical result: {}".format(result["analytical"]))
```

    Analytical result: 150.31276439366033


## Wolfram alpha answers quastion in natural languaje
- What is the average temperature in Bogota in September?




Edición



480b04775e7fac44bb5028ea931a8e3752e75280
ipynb
Jupyter Notebook
SymPyLab/SymPyLab5.ipynb
dgarzona/AlgorithmsUN2020II
ec51783dba40ef930c16040eb1f780c1e53cee97
MIT
SymPyLab/SymPyLab5.ipynb
dgarzona/AlgorithmsUN2020II
ec51783dba40ef930c16040eb1f780c1e53cee97
MIT
SymPyLab/SymPyLab5.ipynb
dgarzona/AlgorithmsUN2020II
ec51783dba40ef930c16040eb1f780c1e53cee97
MIT
Qwen/Qwen-72B
1. YES
2. YES