Proponemos y justificamos un método numérico para calcular la doble integral con límites superiores variables que conduce a la variabilidad de la región de integración. La imposición de variables simples como funciones para los límites superiores proporciona la forma de triángulos de la región de integración y una variable en el límite externo de la integral conduce a un conjunto continuo de triángulos similares. Se superpone una cuadrícula variable en la región de integración. Consideramos tres casos de cambios de la cuadrícula para la división de la región de integración en volúmenes elementales. El primero es solo el tamaño de la cuadrícula impuesta cambia con el cambio de variable del límite superior externo. El segundo caso es el número de elementos de división cambia con el cambio de la variable del límite superior externo. En el tercer caso, el tamaño de la cuadrícula y el número de elementos de división cambian después de fijar su multiplicación. En estos casos, las fórmulas para calcular las integrales dobles se obtienen basadas en la aplicación de cuadraturas en la región interna de integración y realizando la división de triangulación a lo largo del límite variable. El error del método se determina mediante la expansión de la integral doble en la serie de Taylor utilizando el teorema de Barrow. La prueba de eficiencia y confiabilidad de las fórmulas obtenidas del método numérico para tres casos de formas de división de la región de integración se lleva a cabo en ejemplos de la integración doble de funciones suficientemente simples. El análisis de los resultados obtenidos muestra que los errores absolutos y relativos más pequeños se obtienen en el caso de un aumento en el número de elementos de división cuando el aumento de la variable del límite superior externo y el tamaño de la cuadrícula están fijos.