Update README.md

Updated multivariate function and associated figure.

README.md

@@ -11,11 +11,11 @@ McCormick.jl is a component package in the EAGO ecosystem and is reexported by [
 
 Each McCormick object is associated with a
 parameter `T <: RelaxTag` which is either `NS` for nonsmooth relaxations ([Mitsos2009](https://epubs.siam.org/doi/abs/10.1137/080717341), [Scott2011](https://link.springer.com/article/10.1007/s10898-011-9664-7)), `MV` for multivariate relaxations ([Tsoukalas2014](https://link.springer.com/article/10.1007/s10898-014-0176-0), [Najman2017](https://link.springer.com/article/10.1007/s10898-016-0470-0)),
-and `Diff` for differentiable relaxations ([Khan2016](https://link.springer.com/article/10.1007/s10898-016-0440-6), [Khan2018](https://link.springer.com/article/10.1007/s10898-017-0601-2), [Khan2019](https://www.tandfonline.com/doi/abs/10.1080/02331934.2018.1534108)). Conversion between `MV`, `NS`, and `Diff` relax tags are not currently supported. Convex and concave envelopes are used to compute relaxations of univariate functions.
+or `Diff` for differentiable relaxations ([Khan2016](https://link.springer.com/article/10.1007/s10898-016-0440-6), [Khan2018](https://link.springer.com/article/10.1007/s10898-017-0601-2), [Khan2019](https://www.tandfonline.com/doi/abs/10.1080/02331934.2018.1534108)). Conversion between `NS`, `MV`, and `Diff` relax tags are not currently supported. Convex and concave envelopes are used to compute relaxations of univariate functions.
 
 ## **Supported Operators**
 
-In addition to supporting the implicit relaxation routines of ([Stuber 2015](https://www.tandfonline.com/doi/abs/10.1080/10556788.2014.924514?journalCode=goms20)), this package
+In addition to supporting the implicit relaxation routines of [Stuber 2015](https://www.tandfonline.com/doi/abs/10.1080/10556788.2014.924514?journalCode=goms20), this package
 supports the computation of convex/concave relaxations (and associated subgradients) for
 expressions containing the following operations:
 
@@ -45,7 +45,7 @@ In order to bound a function using a McCormick relaxation, you first construct a
 structure that bounds the input variables and then you pass these variables
 to a function.
 
-In the example below, convex/concave relaxations of the function `f(x) = x * (x-5.0) * sin(x)`
+In the example below, convex/concave relaxations of the function `f(x) = x(x-5)sin(x)`
 are calculated at `x = 2` on the interval `[1,4]`.
 
 ```julia
@@ -70,45 +70,51 @@ cv = fMC.cv              # convex relaxation
 cc = fMC.cc              # concave relaxation
 cvgrad = fMC.cv_grad     # subgradient/gradient of convex relaxation
 ccgrad = fMC.cc_grad     # subgradient/gradient of concave relaxation
-Iv = fMC.Intv           # retrieve interval bounds of f(x) on Intv
+Iv = fMC.Intv            # retrieve interval bounds of f(x) on Intv
 ```
 
 Plotting the results, we can easily visualize the convex and concave
 relaxations, interval bounds, and affine bounds constructed using the subgradient
-at the middle of X.
+at the middle of *X*.
 
 ![Figure_1](Figure_1.png)
 
-This can readily be extended to multivariate functions as shown below
+This can readily be extended to multivariate functions, for example, `f(x,y) = (4-2.1x^2+(x^4)/6)x^2+xy+(-4+4y^2)y^2`:
 
 ```julia
+using McCormick
 
-f(x) = max(x[1],x[2])
-
-x = [2.0 1.0]                                    # values of independent variable x
-Intv = [Interval(-4.0,5.0), Interval(-5.0,3.0)]  # define intervals to relax over
-
-# create McCormick object
-xMC = [MC{2,Diff}(x[i], Intv[i], i) for i=1:2)]
-
-fMC = f(xMC)            # relax the function
-
-cv = fMC.cv              # convex relaxation
-cc = fMC.cc              # concave relaxation
-cvgrad = fMC.cv_grad     # subgradient/gradient of convex relaxation
-ccgrad = fMC.cc_grad     # subgradient/gradient of concave relaxation
-Iv = fMC.Intv            # retrieve interval bounds of f(x) on Intv
+# initialize function
+f(x,y) = (4.0 - 2.1*x^2 + (x^4)/6.0)*x^2 + x*y + (-4.0 + 4.0*y^2)*y^2
+
+# intervals for independent variables
+n = 30
+X = Interval{Float64}(-2,0)
+Y = Interval{Float64}(-0.5,0.5)
+xrange = range(X.lo,stop=X.hi,length=n)
+yrange = range(Y.lo,stop=Y.hi,length=n)
+
+# differentiable McCormick relaxation
+for (i,x) in enumerate(xrange)
+    for (j,y) in enumerate(yrange)
+        z = f(x,y)                  # calculate function values
+        xMC = MC{1,Diff}(x,X,1)     # differentiable relaxation for x
+        yMC = MC{1,Diff}(y,Y,2)     # differentiable relaxation for y
+        fMC = f(xMC,yMC)            # relax the function
+        cv = fMC.cv                 # convex relaxation
+        cc = fMC.cc                 # concave relaxation
+    end
+end
 ```
 
-![Figure_3](Figure_3.png)
+![Figure_4](Figure_4.png)
 
 ## Citing McCormick.jl
 
 McCormick.jl is a component of the EAGO.jl ecosystem. Please cite the following paper when using McCormick.jl:
 
 ```
- M. E. Wilhelm & M. D. Stuber (2020) EAGO.jl: easy advanced global optimization in Julia,
- Optimization Methods and Software, DOI: 10.1080/10556788.2020.1786566
+M. E. Wilhelm & M. D. Stuber (2022) EAGO.jl: easy advanced global optimization in Julia, Optimization Methods and Software, 37:2, 425-450, DOI: 10.1080/10556788.2020.1786566
 ```