One possible approach is to decompose each set into a non-overlapping set of axis-aligned rectangles (henceforth boxes), where each edge is shared by at most two rectangles.
Let's say we have a set of two partially overlapping boxes. You split them into non-overlapping sub-boxes (cells), using a grid where each cell is rectangular, but each column of cells has their own width, and each row of cells has their own height:
Here,
$$\begin{array}{lll}
L_1 = x_{min,a} & C_1 = x_{min,b} & R_1 = x_{max,a} \\
L_2 = x_{min,a} & C_2 = x_{min,b} & R_2 = x_{max,b} \\
L_3 = x_{min,b} & C_3 = x_{max,a} & R_3 = x_{max,b} \\
\; & C_4 = x_{max_a} & \; \end{array}$$
$$\begin{array}{lll}
T_1 = y_{min,a} & M_1 = y_{min,b} & B_1 = y_{max,a} \\
T_2 = y_{min,a} & M_2 = y_{min,b} & B_2 = y_{max,b} \\
T_3 = y_{min,b} & M_3 = y_{max,a} & B_3 = y_{max,b} \\
\; & M_4 = y_{max,a} & \; \end{array}$$
Each cell in the grid has width, height, and a color: unused, set 1, set 2, or intersection.
Each vertical cell edge has three pieces of information: its $x$ coordinate, its $height$, and the name of the variable its $x$ coordinate depends on (so that edges defined by a specific variable can be found). Similarly, each horizontal cell edge has an $y$ coordinate, $width$, and the name of the variable its $y$ coordinate depends on.
With the above information, it is trivial to compute the areas. Implement a procedure that sums (width×height) of each cell based on color. Then, the union is the sum of set 1, set 2, and intersect areas; and of course the area of the intersection is just the intersect sum.
More importantly, you can just as easily compute the gradient (the partial derivatives of the area function with respect to each original variable), split by color pairs:
For each horizontal edge of width $w$, examine the color of the cell above, and the color of the cell below. This edge affects the derivative of union with respect to the variable $z$ related to the edge by $dU/dz$, and/or of intersection by $dO/dz$:
- unused, set1: $dU/dz = -w$
- unused, set2: $dU/dz = -w$
- unused, intersection: $dU/dz = -w, \; dO/dz = -w$
- set1, intersection: $dO/dz = -w$
- set2, intersection: $dO/dz = -w$
- intersection, set1: $dO/dz = +w$
- intersection, set2: $dO/dz = +w$
- intersection, unused: $dU/dz = +w, \; dO/dz = +w$
- set1, unused: $dU/dz = +w$
- set2, unused: $dU/dz = +w$
Similarly for the vertical edges.
Some background math to explain this might be in order.
The area of the union $U(x_1, y_1, ..., x_N, y_N)$ is piecewise linear with respect to each coordinate. That is, when one single coordinate (edge of an axis-aligned rectangle) changes by a differential amount $dx_i$, the area of the union changes by amount $(h_{R,i} - h_{L,i})dx_i$, where $h_{R,i}$ is the total length of outside edges defined by this variable on the right side, and $h_{L,i}$ the total length of outside edges defined by this variable on the left side. (This is because increment of the variable on the right edge increments the area, but on the left edge decrements the area.)
This is easiest to understand by looking at a single axis-aligned rectangle, defined by $x_L \le x_R$ and $y_B \le y_T$, where $(x_L , y_B) - (x_R , y_T)$ and $(x_L , y_T) - (x_R , y_B)$ are the two diagonals of the rectangle. Then,
$$A(x_L , y_B , x_R , y_T ) = ( x_R - x_L ) ( y_T - y_B )$$
and the partial derivatives (defining $\nabla A$) are
$$\frac{d\,A}{d\,x_L} = - ( y_T - y_B )$$
$$\frac{d\,A}{d\,y_B} = - ( x_R - x_L )$$
$$\frac{d\,A}{d\,x_R} = ( y_T - y_B )$$
$$\frac{d\,A}{d\,y_T} = ( x_R - x_L )$$
When we decompose the set of overlapping boxes to non-overlapping boxes in a grid, the above applies to each box, if and only if the other box sharing the same edge is of a different type: if it were of the same type, moving the single edge between the two boxes by an infinitesimal amount would not have any effect on their total area at all.
This is not a particularly hard programming problem, because there are many different ways to implement the data structures needed to solve it. (Although it does mean that finding a near-optimal solution is hard, simply because there are so many different ways of implementing this, and only actual practical testing would show which ones are most effective.)
I would suggest implementing it and testing it separately first, perhaps having the test program output an SVG image of the resulting decomposed non-overlapping set, with outside horizontal edge widths and vertical edge heights and the variable names their coordinates depend on written on the outside of the image, for area and gradient verification by hard.
If this approach feels intractable, it should be noted that for $N$ rectangles ($4N$ independent variables), you could simply calculate the derivatives using
$$\frac{d\,A(c_1, c_2, ..., c_{4N})}{d\,c_i} = \frac{A(..., c_i + \Delta_i, ...) - A(..., c_i - \Delta_i, ...)}{\Delta_i}$$
(involving $8N$ area calculations, so linear complexity with respect to $N$), where $\Delta_i$ is a small perturbation in the coordinate $c_i$. In a very real sense, if $\Delta_i$ approximates the typical change in one iteration step (if done along variable $c_i$), this should give more realistic feedback! You see, since the area functions are piecewise linear, and there are up to $2N-1$ pieces along each axis, the derivative, or infinitesimal or instantaneous change along that axis, may not be truly informative!
For example, consider a case where the two sets have a long thin hole somewhere. Because the hole has a long edge, the derivative with respect to the variables defining the long edges of the holes are large, and that makes those variables "more important" when looking at the derivatives only. In reality, the area of the hole may be minuscule compared to the totality, making those variables not at all important in reality. If $\Delta_i$ for the variables is large enough to span the hole, the estimate obtained with it better reflect the actual change in area if the hole were covered.
Best Answer
It is not always possible to attain $0.359$.
Assuming two congruent rectangles of base $1$ and height $h$, note that if the rectangles coincide, or are concentric at a $90^o$ orientation, we can have two disjoint rectangles each equal to $\frac{h}{2}$.
But if the rectangles are concentric at $45^o$, it appears that greatest equal disjoint rectangles can be less than $.359$ of the containing rectangles for$$1<h<2$$To verify this we must generalize somewhat the equation offered by @Oscar Lanzi.
In the figure below, concentric congruent rectangle $FGHI$ is rotated $45^o$ from rectangle $BCDE$ about center $A$. If rectangle $BJ$ within $BCDE$ is equal to rectangle $MNHO$ within $FGHI$, then denoting $LM$ as $x$ and $LH$ as $l$, we can write (cf. @Oscar Lanzi)$$hx=\frac{1}{2}(l-x)^2$$ Note we are assuming $MH$ is a square, since we are seeking whether some $x$ yields a pair of greatest equal, non-overlapping rectangles $BJ$, and $MH$ with vertex at $H$ and opposite vertex $M$ on the right side of $BJ$, each of which is less than $.359$ of the rectangle that contains it. For each lesser $BJ$ there is of course some equal disjoint rectangle $MH$, but we seek a lower limit for the greatest equal rectangles. And Euclid Elements, VI, 27 shows that, within isosceles right triangle $KHQ$, square $MNHO$ is the greatest rectangle bounded by the equal sides and the hypotenuse.
But since the diagonal $FH$ of $FGHI$ lies along $LM$ only when the rectangles are squares, for Lanzi's $l=\frac{1+\sqrt2}{2}$ we must write the more general expression$$l=\frac{1+cos\angle AHP\sqrt{1+h^2}}{2}$$where $AP$ is drawn parallel to $BE$.
For while $LP=\frac{1}{2}$$$PH\ne\frac{1}{2}FH=AH=\frac{\sqrt{1+h^2}}{2}$$but rather$$PH=AHcos\angle AHP$$
Hence we have$$hx=\frac{1}{2}[(\frac{1+cos\angle AHP\sqrt{1+h^2}}{2}-x]^2$$which yields finally the quadratic$$x^2+(-2h-1-cos\angle AHP\sqrt{1+h^2})x+\frac{1+2cos\angle AHP\sqrt{1+h^2}+cos^2\angle AHP(1+h^2)}{4}=0$$whose one solution (since $x<1$) is$$x=\frac{2h+1+cos\angle AHP\sqrt {1+h^2}}{2}-\sqrt {h^2+h+hcos\angle AHP\sqrt{1+h^2}}$$Example 1. For $h=\sqrt2$ we have $FH=\sqrt3$, giving $AH=\frac{\sqrt3}{2}$, whence $\angle FHG=arccos\frac{\sqrt2}{\sqrt3}=35.26^o$, and $\angle AHP=45-35.26=9.74^o$, whose $cos=.986$.
Substitution then gives$$x=2.768-2.414=.354$$ And since areas of rectangles under the same height are proportional to their bases$$BJ=.354\times BCDE$$and we have two greatest equal disjoint rectangles whose areas are less than .359 of the congruent rectangles containing them.
Example 2. Again, when $h=\sqrt3$ and consequently $\angle AHP=15^o$, we find by the same method that$$x=3.198-2.842=.356$$Example 3. But for $h=2$, when $\angle AHP=18.47^o$$$x=3.561-3.201=.360$$And for $h>2$ the excess over $.359$ increases.
But note that $h$ soon reaches a limit, since for $h\approx2.652$ we get $x=.377$, making$$hx=BJ=MNHO=1$$ and this is the greatest value $MNHO$ can have and still lie within rectangle $FGHI$.
The case of rectangles generally remains to be considered, but the foregoing seems to show that at least within certain congruent rectangles there can be greatest equal disjoint rectangles whose areas are less than $.359$ of the rectangles that contain them.