Visualization of the convex integration solutions to the Monge-Amp\`ere equation

In this article, we implement the algorithm based on the convex integration result proved in [Lewicka-Pakzad Analysis and PDE (2017)] and obtain visualizations of the first iterations of the Nash-Kuiper scheme, approx- imating the anomalous solutions to the Monge-Amp\`ere equation in two dimensions.


Introduction
The purpose of this paper is to implement the technique of convex integration to obtain and analyze the visualizations of the Hölder regular C 1,α solutions to the Monge-Ampère equation: We rely on the results of the paper [13], where the Nash-Kuiper iteration was shown to be adaptable to prove flexibility of the very weak solutions to (1.1) in the regularity regime 0 < α < 1 7 . More precisely, the following holds: Theorem 1.1. [13] Let f ∈ L 7/6 (Ω) on an open, bounded, simply connected domain Ω ⊂ R 2 . Fix an exponent α ∈ (0, 1 7 ). Then the set of functions v ∈ C 1,α (Ω) that satisfy (1.1) in the sense of distributions, is dense in the space C 0 (Ω). Namely, for every v 0 ∈ C 0 (Ω) there exists a sequence v n ∈ C 1,α (Ω) converging uniformly to v 0 and solving (1.1). When f ∈ L p (Ω) and p ∈ (1, 7 6 ), the same result is true for any exponent α ∈ (0, 1 − 1 p ). We note that by a straightforward approximation argument, if v ∈ W 2,2 then Det ∇ 2 v = det ∇ 2 v a.e. in Ω, where ∇ 2 v ∈ L 2 (Ω) stands for the standard Hessian matrix field of v. It has been shown [12,18] that in this case, condition f ≥ c > 0 in Ω (respectively, f ≡ 0 in Ω) implies that v must be C 1 and globally convex or concave (respectively, developable). Likewise [13,14], the flexibility result of Theorem 1.1 has its rigidity counterpart: if v ∈ C 1,α (Ω) with α > 2 3 solves (1.1) for the right hand side that is positive Dini continuous (respectively, f ≡ 0 in Ω) then v is a convex Alexandrov solution to det ∇ 2 v = f (respectively, v is developable).
For a better understanding of the above statements and their relation to nonlinear elasticity of thin films, we refer to [13]. We now recall the connection between solutions to (1.1) and the problem of isometric immersion of Riemannian metrics. Consider the one-parameter family of metrics h = Id 2 + 2 2 A on Ω, where A : Ω → R 2×2 sym is a perturbation tensor that satisfies: (1.2) − curl curl A = f, for example one can take: A = −∆ −1 (f )Id 2 . Consider also the pull-back metrics (∇φ ) T ∇φ of the family of deformations φ = id + ve 3 + 2 w : Ω → R 3 . Computing the Gauss curvature of h : and the curvature of the pull-back metrics: κ((∇φ ε ) T ∇φ ε ) = − 2 curl curl 1 2 ∇v ⊗∇v +sym ∇w + o( 2 ), we see that (1.1) can be interpreted as requesting the out-of-plane displacement v to be amenable to matching by a higher order in-plane displacement w, to achieve the prescribed Gaussian curvature f , at its highest order term. Indeed, on a simply connected domain the kernel of the differential operator "curl curl" consists of the fields of the form sym ∇w. Thus, equation (1.1) is equivalent to: where A satisfies (1.2) and w : Ω → R 2 is an auxiliary variable (a Lagrange multiplier). In the same spirit, equation (1.3) can be interpreted as the "infinitesimal isometry" constraint, namely: the family h should coincide, up to terms of order 2 , with the pull-back metrics (∇φ ) T ∇φ defined above. In order to construct C 1,α solutions of (1.3) with v approximating a given continuous v 0 (by density, it suffices to work with v 0 smooth), it has been shown in [13] that a Nash-Kuiper convex integration method can be employed. This method [7,11] modifies an initial "short infinitesimal isometry" (v 0 , w 0 ), i.e. a couple for whom (1.3) is satisfied with an inequality (hence a subsolution) rather than the equality, towards an exact solution, in successive small steps. Technical similarities with the isometric immersion construction [16,17,11,4] are based on the presence of the quadratic terms (∇φ ) T ∇φ and ∇v ⊗ ∇v in both PDEs. For parallel application of the convex integration technique to constructing the energy-dissipative solutions of the Euler equations and to the recent resolution of the Onsager's conjecture, we refer to [5,6,3,8].
As we have mentioned, Theorem 1.1 follows from the main result in [13] below: Theorem 1.2. [13] Let Ω ⊂ R 2 be an open and bounded domain. Let v 0 ∈ C 1 (Ω), w 0 ∈ C 1 (Ω, R 2 ) and A ∈ C 0,β (Ω, R 2×2 sym ), for some β ∈ (0, 1), be such that: Then, for every exponent α ∈ 0, min 1 7 , β 2 , there exist sequences: v n ∈ C 1,α (Ω) which converges uniformly to v 0 , and w n ∈ C 1,α (Ω, R 2 ) which converges uniformly to w 0 , satisfying: In this paper, we implement the proof of Theorem 1.2 (for f and v 0 as in Theorem 1.1 and specially prepared A and w 0 ) into an explicit algorithm, leading to the visualizations of the anomalous solutions to the Monge-Ampère equation (1.1). These are solutions failing to possess the indicated geometrical structural properties valid in the rigid regime v ∈ W 2,2 or v ∈ C 1,2/3+ : • nonconvex solutions, approximating a saddle v 0 (x, y) = x 2 − y 2 , when f ≡ 1, • nondevelopable solutions, approximating a convex paraboloid v 0 (x, y) = x 2 + y 2 , when f ≡ 0, • solutions uniformly close to v 0 (x, y) = x 2 + y 2 , when f ≡ −1. We remark that numerical implementation of convex integration construction for the isometric immersion problem has been reported in [1,2], leading among others to the remarkable images of flat torus embedded in R 3 .
The paper is organized as follows. In sections 2, 4 and 5 we revisit the proof of Theorem 1.2 and explicitly trace all the constants appearing in the arguments in [13] for an optimal performance of the numerical algorithm. Without loss of generality, we will assume that 0 ∈ Ω. In sections 3 and 6 we implement the indicated technique and visualise the first few approximations of anomalous solutions to (1.1) in case of f being positive, negative and equivalently equal 0, respectively. The implementation is done in MATLAB, where for the detailed Hölder continuous approximations in section 6 we use of the Python package mpmath [15], allowing the user to define floating point arithmetics up to an arbitrary precision. 1406730 and DMS-1613153.
2. The C 1 convergence and construction of corrugations We start with a weaker version of Theorem 1.2. Define three unit vectors: , we have the following. For any > 0 there exist v ∈ C 1 (Ω) and w ∈ C 1 (Ω, R 2 ) that satisfy: Moreover, if with some smooth positive functions φ k ∈ C ∞ (Ω) there holds: then it is possible to have the field w in (2.1) additionally satisfy: w − w 0 0 ≤ . The vector field (v, w) is obtained as a C 1 (Ω)-limit of smooth fields v k ∈ C ∞ (Ω) and w k ∈ C(Ω, R 2 ) such that v k − v 0 0 ≤ and w k − w 0 0 ≤ for all k ≥ 1.
As noted below, condition (2.2) implies positive definiteness of the left hand side: We now recall the convex integration construction which allows to decrease such positive definite defect, in each of the η k ⊗η k directions, through oscillatory perturbations in v and w. Let η be a unit vector and φ a smooth positive function onΩ. Given two fields v and w, we want to adjust them toṽ andw so that: Firstly, defineṽ .
= v + 1 λ f (x, λx · η) by adding oscillations of amplitude 1 λ and frequency λ, with λ → ∞ in the ultimate limiting process. Secondly, decompose perturbation of w in the two directions ∇v and η, and writew .
, and h(x, t) should be bounded and 1-periodic in the variable t. The expansion: The actual choice of f, g and h is given in Lemma 2.2.
In the subsequent proof of Theorem 2.1, we will aim at keeping the frequencies λ = λ k minimal, facilitating a good numerical implementation. Three modifications with respect to proofs in [13] will be done in order to improve the estimates on λ k : (i) all the global inequalities will be specified to pointwise estimates; (ii) we will allow the error threshold parameter δ to depend on x rather than be constant; (iii) we will set distinct rates of convergence for the errors D k and the approximations v k in Proposition 2.5.
Proof. The three indicated rank-one matrices are linearly independent by a straightforward calculation. The formula in (i) follows directly as well and implies, in view of the Cauchy-Schwartz inequality: |φ 1 | ≤ √ 17 4 |B|. In the same manner, we have: In the setting of (iv), since b kk ≤ |B| for k = 1, 2 and b 12 ≤ √ 2 2 |B|, we get:φ with some φ k = φ k (x) and a constant d > 0. Let ξ > 0 be such that: Fix > 0; then there existṽ ∈ C ∞ (Ω),w ∈ C ∞ (Ω, R 2 ) and a constantd > 0 such that: where C is a universal constant.
. Define smooth, positive functions δ, a k :Ω → R: Clearly, (2.13) guarantees the bound: , the successive corrections v k and w k are now constructed by applying Proposition 2.3 to v = v k−1 , w = w k−1 , a = a k , η = η k and an appropriate λ = λ k ≥ 1 determined below in (2.21) and (2.23). Observe that: where (2.7) yields the following pointwise bound on the matrix-valued fields {B k } 3 k=1 :

2.
To prove positivity of the decomposition in (2.14), we set: and use Lemma 2.4 to: 3 Namely, by Lemma 2.4 (ii) it follows that for k = 1 . . . 3 we have: where the second inequality above is valid when: Note that the first estimate in (2.15) holds then as well, again in view of Lemma 2.4 (ii): For the validity of (2.20), we choose {λ k } 3 k=1 so that: (2.21) 3. Observe that, by Lemma 2.4 and the definition of a k , we get: By (2.8), there follows the second inequality in (2.15): if only we assume the following condition: Note that the third inequality in (2.21) and the bound on φ k in terms of |D| yield: Thus, by (2.9) and (2.22) we conclude the first bound in (2.16): Using the second inequality in (2.8) together with (2.23) and (2.22), we obtain the third bound in (2.15): whereas (2.9) and (2.24) are used for the final estimate in (2.16): Above, the term has been bounded by 1 2 |D| 0 in view of (2.21).
Proof of Theorem 2.1.
1. Assume first that (2.2) holds. We will construct a sequence of smooth approximations which converge in C 1 to the required solution in (2.1). Starting with v 0 , w 0 , we define recursively v k+1 and w k+1 by applying Proposition 2.5 to v k ∈ C ∞ (Ω) and w k ∈ C ∞ (Ω, R 2 ), where we denote the corresponding defect D k . = A − ( 1 2 ∇v k ⊗ ∇v k + sym∇w k ), and request that the construction parameters k , ξ k > 0 satisfy: for an appropriately small constant β ≤ 1. By (2.14), each D k can be decomposed in the basis with strictly positive coefficients. Further, (2.15) implies: while by (2.16) and (2.15) we obtain: Consequently, the sequence {v k } ∞ k=0 is Cauchy in C 1 (Ω) and thus it converges to some v ∈ C 1 (Ω). By (2.26), the first statement of (2.1) follows. In particular, the norms ∇v k 0 are uniformly bounded (by C( ∇v 0 0 + D 0 0 + 1), which allows to compute: where C is a universal constant. Hence, {w k } ∞ k=0 is Cauchy and converges in C 1 (Ω, R 2 ) to some w ∈ C 1 (Ω, R 2 ), satisfying: This proves the second statement in (2.1), by (2.15) and (2.25). Also, it is clear that taking

2.
In the absence of the positivity of the decomposition , the decomposition of the modified defect: Applying the first part of the proof to v 0 ,w 0 andD 0 yields v ∈ C 1 (Ω) and w ∈ C 1 (Ω, R 2 ) with the desired properties.

A numerical implementation of Theorem 2.1
The purpose of this section is to obtain images of the first few approximation steps in the convex integration construction of the solutions to the Monge-Ampère equation (1.1). We will consider two case scenarios: (i) the right hand side f (x, y) = 1 and the subsolution v 0 (x, y) = x 2 − y 2 is non-convex; (ii) the right hand side f (x, y) = −1 and the subsolution v 0 (x, y) = x 2 + y 2 is strictly convex on the domain Ω = (−0.5, 0.5) × (−0.5, 0.5). We set the threshold = 0.1 in seeking a solution v with v − v 0 0 < 0.1 as in Theorem 2.1.
We were able to exhibit three consecutive corrugations. Since the frequencies {λ k } quickly become very large, fine meshes needed to be used, requiring a great computing power. On the other hand, a priori estimates (2.21) and (2.23) could be efficiently validated; indeed these estimates helped to experimentally reduce the values of each λ k .
For the first approximation we sampled on a square grid with the initial step size 0.001, followed by the step size of the order h ∼ 0.1 λ k . With this choice of h we were so far able to perform two steps of convex integration in both examples below. The derivatives on the mesh are: We assigned the values of the initial auxiliary variable w 0 to make the initial defect At each step, the defect D was calculated explicitly and then decomposed into D = 3 k=1 φ k η k ⊗ η k using the formulas from Lemma 2.4. We set δ = 0.5 and followed the indicated construction. By (2.16), our numerical implementation reduced the norm of defect to 3 4 of the previous defect norm every three convex integration steps.
Define the coefficients:ã k = φ k 2 and the auxiliary matricesB k by: while in order to getφ k ≥d = d 4 = 0.1, we proceed as in (2.19) and introduce the condition: Conditions (3.1) and (3.2) were used to determine the frequencies λ k , where in estimating the three terms in eachB k according to (2.18), we bound |∇ 2 v k−1 | using (2.10) with the previous choice of λ k−1 . Once the three choices of λ k were made, we performed a verification that the new v and w obtained with the described above modification step remained indeed within = 0.1 error from the original fields.
The values of λ 1 thus obtained were small enough to allow for the numerical executing the first step of convex integration. We defined the fields v 1 and w 1 as in Proposition 2.3, calculatedB 1 , its decomposition intoB 1 = 3 k=1 φ k η k ⊗ η k , and verified the two conditions (3.1) and (3.2). This operation has been repeated until the smallest possible λ 1 was found. We then used this value and the resulting a priori estimates to obtain the second frequency λ 2 , small enough to apply the second step of convex integration. Thus we obtained the fields v 2 and w 2 , and, by the same procedure described above, the third optimal frequency λ 3 together with v 3 and w 3 . The final step size h needed to be reduced to 0.0001 in order to allow 10 mesh points for each corrugation.
Example 3.1. We approximate v 0 (x, y) = x 2 − y 2 with a solution v to: Take w 0 (x, y) = (xy 2 , x 2 y) and A(x, y) = 5 − x 2 +y 2 4 Id 2 , so that: −curl curl A(x, y) = 1. Then: and the corresponding defect is diagonal and positive definite: The function v 0 and its two subsequent corrugations are shown in Figure 3.1. In Figure 3.2 we show a detailed picture of the second corrugation; the red portion of the graph indicates the area on which we applied the third corrugation shown in Figure 3   Example 3.2. We approximate v 0 (x, y) = x 2 + y 2 with a solution v to: In this example we take w 0 (x, y) = (−xy 2 , −x 2 y) and A(x, y) = 5 + x 2 +y 2 4 Id 2 , satisfying −curl curl A(x, y) = −1 and resulting in the diagonal, positive definite defect: We plot three images starting from v 0 and subsequently adding the first and second corrugations in Figure 3  give an upper estimate of the uniform distance of v obtained through three steps of convex integration from the initial subsolution v 0 . The value ( B 1 0 + B 2 0 )/ D 0 which does not take the third corrugation into account, needs to be below 2 · 1 12 = 1 6 . The contribution of the last corrugation is guaranteed to be less than 1 12 D 0 through the a priori estimates. Lastly, min φ k are the minima of each of the coefficients inΩ, in the defect computed after two corrugations. The a priori estimates once again guarantee that the third step will not make these less than the required error value 0.1.     In this and the next section we give a proof of Theorem 1.1 using a constructive, numerically implementable algorithm. The details of the implementation and the resulting visualizations will be presented in section 6. The construction follows the proof in [13], we however make some important modifications and compute all the relevant constants explicitly. We begin with a choice of a standard mollifier and some preliminary estimates.
Let φ ∈ C ∞ c (B(0, 1)) be the following radially symmetric function: where A =´B 1 (0) exp(− 1 1−|x| 2 ) dx = π · (1/e + Ei(−1)) so that´R 2 ϕ = 1. The constant A, given in terms of the exponential integral Ei, may be approximated to any degree of precision. We will use A ∈ [0.46, 0.47] in the estimates below, but when implementing the algorithm numerically we will evaluate A more precisely.
The next result is a modification of Proposition 2.3 so we omit its proof.
Proof. The proof proceeds in three parts: mollification of v, w, A to control higher derivatives, modification of w to ensure the positive decomposition of the defect in the basis {η k ⊗ η k } 3 k=1 , and application of three consecutive steps of convex integration to reduce the defect.

Mollification.
For all l < r, define the following functions onΩ through mollification with the standard kernel ϕ as in (4 .1): Denote: D . = A − 1 2 ∇v ⊗ ∇v + sym∇w . We use Lemma 4.1 and assumption (4.11) to obtain the uniform error bounds below, where the relevant norms of quantities v, w, A, D are taken on Ω, while other norms are taken on the supersetΩ r : (4.16) In view of (4.2), (4.6) and (4.11), the last bound above is specified to: (4.17) We also get, by (4.3):

Modification and decomposition.
To ensure that the deficit may be decomposed with positive coefficients, we define: By (4.16) there follow the bounds: 4 , √ 2 + 9 5 and therefore: ∇ mD = ∇ m D for all m ≥ 1. Further, this construction guarantees that in the decompositionD = 3 k=1 φ k η k ⊗ η k on Ω, in view of Lemma 2.4 (iv) we get: φ k ≥ ( D C 0 (Ωr) + D 0 )/2. We now find the bounds on the first norms of the smooth positive functions a k . = √ φ k . We begin by noting that min x∈Ω a k (x) ≥ . In view of (4.20), (4.17) and since ∇a k = ∇φ k 2a k and ∇ 2 a k = ∇ 2 φ k 2a k − ∇a k ⊗∇a k a k , Lemma 2.4 yields: (4.21) 3. Iteration of convex integration. We set v 0 .
= v, w 0 . =w restricted toΩ and then define recursively v k ∈ C 3 (Ω), w k ∈ C 2 (Ω, R 2 ) for k = 1, 2, 3 by applying Proposition 4.2. We apply it to v k−1 , w k−1 , with a k from the decomposition ofD in the basis given by {η k } 3 k=1 . Lastly, we set the parameters: and the non-decreasing triple {δ k } 3 k=1 with the initial choice: The construction is complete by settingṽ .

A proof of Theorem 1.1
We start by showing the main approximation result needed in Theorem 1.1, that is a version of the result in Theorem 1.2. It consists of iterating the "stages" construction, with the sole restrictive assumption on the smallness of the initial deficit. Given three functions: v ∈ C 2 (Ω r ), w ∈ C 2 (Ω r , R 2 ) and A ∈ C 0,β (Ω r , R 2×2 sym ) of Hölder regularity β ∈ (0, 1), assume that: Then, for every exponent α ∈ (0, min 1 7 , β 2 ) one can findv ∈ C 1,α (Ω),w ∈ C 1,α (Ω, R 2 ) such that: Proof. 1. Define a decreasing sequence of open domains {Ω k } ∞ k=0 by setting: , resulting in the nonzero deficit: we will construct the functions v k+1 ∈ C 2 (Ω k+1 ), w k+1 ∈ C 2 (Ω k+1 , R 2 ) by applying Proposition 4.3 to the sets Ω k+1 ⊂ Ω k , the matrix field A |Ω k and the parameters σ k , M k chosen according to the procedure indicated below. First, choose the exponent s ∈ (0, 1) to satisfy: Existence of such s in guaranteed by α ∈ (0, min 1 7 , β 2 ). Second, set the constant: and let {σ k } ∞ k=0 be an increasing sequence of positive numbers, converging to some σ max and satisfying: We note that it is enough to take σ k = σ max for all k, but having σ k as small as possible is advantageous for the numerical calculations. Third, the constants {M k } ∞ k=0 are defined by: for a large constant N ≥ 1 below (that is evaluated numerically in the following implementation). In particular, there holds: 2. We now inductively prove that: For k = 1, Proposition 4.3 gives: The second term in the right hand side above is majored by 1 2σ s 0 by (5.3), whereas the first term is also bounded by: (5.4). For k ≥ 1, we similarly observe that: where the bound on the second term follows as in the case k = 0, while to estimate the first term: we used the inductive assumption, (5.4) and (5.1), implying that 3β + s( β 2 − 1) > 0. This ends the proof of (5.6).
Observe also that by (5.6), (5.3) and the assumed bound on δ 0 we get, for all k ≥ 0: 3. Clearly, when D k = 0, we stop the recursive process and set (v,w) . = (v k , w k ) with the claimed error estimates established in the next step below. We now proceed by validating the assumptions of Proposition 4.3 in case D k = 0.

4.
We now show that the sequences {v k } ∞ k=0 , {w k } ∞ k=0 converge in C 1,α (Ω) to somev,w. By (5.6) this will imply that A − 1 2 ∇v ⊗ ∇v + sym∇w = lim k→∞ D k = 0 together with the desired estimates on v − v 0 and w − w 0 . Indeed: by Proposition 4.3, (5.6), (5.4). Automatically, the bound above implies that {v k } ∞ k=0 is Cauchy in C 0 (Ω). The similar statements for {w k } ∞ k=0 follow by (5.7) in: Finally, we estimate the C 1,α norms by interpolating from C 1 to C 2 norms in: Above, the constant C depends, in its first appearance, only on the curvature of a smooth superset ofΩ contained in Ω r . Eventually, C depends on various numerical values including M 0 , α and σ max but it is independent of k. We see that the last condition in (5.3) implies that the constant of the geometric progression in the right hand sides is less than 1 for large k, that is when σ k is close to σ max . This establishes that {v k } ∞ k=0 , {w k } ∞ k=0 are Cauchy sequences in C 1,α (Ω) and completes the proof of Proposition 5.1.

Proof of Theorem 1.2.
Let f, v 0 be as in Theorem 1.1. Fix a small parameter > 0. We first extend v 0 to a continuous function on R 2 and approximate this extension withv 0 ∈ C ∞ (R 2 ) satisfying: where c > 0 is a constant and λ ∈ C 0,β (V ) solves: −∆λ = χΩf in an open superset V ofΩ, with λ = 0 on ∂V . Then, as discussed in the introduction, A = 1 2 ∇v 0 ⊗ ∇v 0 + sym∇w 0 holds inΩ r for some w 0 ∈ C 1 (Ω r , R 2 ) on Ω r = Ω + B r (0) for some r > 0. Note that if c is large enough to guarantee: then in view of Lemma 2.4 it follows that: Clearly, the field w 0 may be exchanged with a smooth fieldw 0 ∈ C ∞ (Ω r , R 2 ) so that the positive coefficient decomposition above is still valid. Applying now Theorem 2.1, we get v ∈ C ∞ (Ω r ) and w ∈ C ∞ (Ω r , R 2 ) such that: v −v 0 0 ≤ and A − 1 2 ∇v ⊗ ∇v + sym∇w 0 ≤ min , r, (5.4)10 −16 .
In virtue of Proposition 5.1, we further findv ∈ C 1,α (Ω) andw ∈ C 1,α (Ω, R 2 ) solving: 1 2 ∇v ⊗ ∇v + sym∇w = A inΩ, with the approximation error: Concluding, the functionv is a Hölder regular solution to (1.1), approximating the given v 0 with the arbitrarily small error: as claimed in the result. 6. Numerical implementation of the C 1,α convergence scheme In section 3 we implemented the approximation construction of Theorem 2.1 and Proposition 2.5; we now implement the Hölder approximation specified in Proposition 4.3 and Proposition 5.1. Visualizations relevant to the C 1,α case are much harder to obtain. The main difficulties come from the smallness of various parameter values and and the fact that the ratio σ between the frequencies of subsequent corrugations is by necessity very large. This first obstacle was overcome through the use of the Python package mpmath [15], allowing the user to define floating point arithmetics up to an arbitrary precision. The second obstacle of the necessity of large σ is insurmountable in most cases. Nonetheless, in the example of the degenerate Monge-Ampére equation det ∇ 2 v = 0, small values of λ were found to produce a reduction in the deficit which justified the visualization obtained below.
The approach taken was significantly different from the case of C 1 convergence. In each of the examples the problem was solved explicitly through symbolic calculations, defining symbolic variables and parameters in MATLAB (the only calculation which could not be done symbolically was the mollification step). The resulting text files, containing the outcomes were then modified through a script to make the syntax compatible with the mpmath package. We subsequently sampled these functions at 1000 randomly selected points on the square domain Ω = (−1, 1) × (−1, 1), and the maximum value attained by each function was recorded. The main purpose of these calculations was to study how the results of convex integration varied by changing the value of σ. We argue that the maximum sample value, while not being the exact supremum of these functions on the domain, is a good approximation of the order of magnitude. In fact, when the process was repeated over different sets of 1000 sample points, the variation in the values obtained was negligible when considering the order of magnitude, as illustrated in the table below. In this table we show the results of running our code 10 times on the same example with the same parameters, where the only change was the choice of the sample points. In each test we evaluated the maximum defect recorded, the maximum value of |v| and the maximum gradient of w. As can be seen, the order of magnitude of all these functions remains the same. Similar results were obtained with different examples, functions and parameters, but the details are omitted as they do not add particular insight.
In each of the examples a value of λ 1 = 10 19 was chosen and σ was increased exponentially from 10 1 to 10 18 , covering the orders of magnitude between the theoretical minimum of σ and σ max . We then applied the modification of w and the three steps of convex integration in Proposition 4.3 and evaluated the maximum of the new defect, the norm of v together with is gradient and Hessian. The numerical results obtained are summarized in the tables in Appendix A. [19] V.Šverák, On regularity for the Monge-Ampère equation without convexity assumptions, preprint, Heriot-Watt University (1991).