# Solve system of equations with integration over regions

I want to solve the following problem in Mathematica:

Assume that $$\lambda_2=-\frac{1}{2}+ r+\frac{1}{1+r}- r^2\log\frac{1+r}{r}$$, $$v\in[0,1]$$ and $$\theta\in[0,1]$$. For pre-specified values of $$k,r>0$$, solve for $$\{s,\lambda_1\}$$ through the following two equations $$\int_{\left\{\stackrel{2\theta v s+\lambda_2(1-\theta)r\geq\lambda_1}{\theta(v+r)\geq r}\right\}}(1-\theta)r\,dvd\theta+ \int_{\left\{\stackrel{2\theta v s+\lambda_2v\geq\lambda_1}{\theta(v+r)\geq r}\right\}}v\theta\,dvd\theta=s$$ $$\int_{\left\{\stackrel{2\theta v s+\lambda_2(1-\theta)r\geq\lambda_1}{\theta(v+r)\geq r}\right\}}\,dvd\theta+ \int_{\left\{\stackrel{2\theta v s+\lambda_2v\geq\lambda_1}{\theta(v+r)\geq r}\right\}}\,dvd\theta=k$$

For the first step, I tried to calculate the integrals using Integrate and ImplicitRegion, but got stuck here because Mathematica refuse to evaluate the expressions. Here are my codes:

λ2 = -(1/2) + r + 1/(1 + r) - r^2 Log[1 + 1/r];
k = 0.5;
r = 0.5;
R1 = ImplicitRegion[θ (v + r) > r && 2 θ v s + (1 - θ) r > λ1, {{θ, 0, 1}, {v, 0, 1}}];
R2 = ImplicitRegion[θ (v + r) < r && 2 θ v s + λ2 v > λ1, {{θ, 0, 1}, {v, 0, 1}}];
Integrate[(1 - θ) r, {θ, v} ∈ R1] + Integrate[θ v, {θ, v} ∈ R2]
Integrate[1, {θ, v} ∈ R1] + Integrate[1, {θ, v} ∈ R2]


Ultimately, I would like to plug back the solved $$\{s,\lambda_1\}$$, and plot the regions of R1 and R2. But how can I fix the codes to solve the equations?

• Mathematica may not be able to solve this problem. However, to improve its chances of doing so, can you provide any bounds on s and λ1? Jul 31, 2020 at 3:12
• Also, supposing that Integrate can be performed, how are s and λ1 to be determined from the result? Jul 31, 2020 at 3:51
• GIven the specified values of k and r, the only unknowns would be s and lambda1, and we have two equations. So, hopefully, that would do the trick. I do anticipate that s is between 0 and 1, and lambda1 is positive. Jul 31, 2020 at 4:48

This interesting problem can be solved numerically by computing InterpolationFunctions for the two sums of integrals in the last two lines of code in the Question.

λ2 = -(1/2) + r + 1/(1 + r) - r^2 Log[1 + 1/r];
k = 0.5;
r = 0.5;

t = Flatten[Table[
R1 = ImplicitRegion[θ (v + r) > r && 2 θ v s + (1 - θ) r > λ1, {{θ, 0, 1}, {v, 0, 1}}];
R2 = ImplicitRegion[θ (v + r) < r && 2 θ v s + λ2 v > λ1, {{θ, 0, 1}, {v, 0, 1}}];
{{s, λ1},
Integrate[(1 - θ) r, {θ, v} ∈ R1] + Integrate[θ v, {θ, v} ∈ R2],
Integrate[1, {θ, v} ∈ R1] + Integrate[1, {θ, v} ∈ R2]},
{s, 0, 1, .05}, {λ1, 0, 1, .05}], 1];


(I assume here that s and λ1 both lie between zero and one. If not, the limits on Table must be adjusted accordingly.)

f1 = Interpolation[Delete /@ t];
Plot3D[f1[s, λ1], {s, 0, 1}, {λ1, 0, 1}, AxesLabel -> {s, λ1}, ImageSize -> Large,
LabelStyle -> {15, Bold, Black}] f2 = Interpolation[Delete /@ t];
Plot3D[f2[s, λ1], {s, 0, 1}, {λ1, 0, 1}, AxesLabel -> {s, λ1}, ImageSize -> Large,
LabelStyle -> {15, Bold, Black}] In the absence of equations to determine s and λ1 in terms of f1 and f2, use {f1[s, λ1] == .03, f2[s, λ1] == .135} for illustrative purposes. Then, the resulting values of s and λ1 are

sol = FindRoot[{f1[s, λ1] == .03, f2[s, λ1] == .135}, {{s, .25}, {λ1, .35}}]
(* {s -> 0.153035, λ1 -> 0.362465} *)


Note that FindRoot requires fairly good initial guesses here to avoid attempting to search outside the domains of f1 and f2. Plots of the corresponding regions are

RegionPlot[(θ (v + r) > r && 2 θ v s + (1 - θ) r > λ1) /. sol, {θ, 0, 1}, {v, 0, 1},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}, PlotPoints -> 60] RegionPlot[(θ (v + r) < r && 2 θ v s + λ2 v > λ1) /. sol, {θ, 0, 1}, {v, 0, 1},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}, PlotPoints -> 60] • Very tricky answer! Where did you find Delete /@ t to remove the third element, never seen it. Thanks. Jul 31, 2020 at 7:16
• @UlrichNeumann I had planned to use Delete[#, 2]& /@ t to remove the second column, then noticed the operator form of Delete in the documentation and decided to use it instead. I used Delete instead of Most so that the code for f1 and f2 would look similar. Several list manipulation functions, for instance Cases, also have operator forms. Best wishes. Jul 31, 2020 at 13:47
• Thanks it' s very helpful! Jul 31, 2020 at 13:59
• Thank you, this is super helpful! I have a syntax question, why do you need to Flatten the table? Jul 31, 2020 at 16:01
• @user391830 Table with two indices produces a nested List. Interpolation expects a flat List. See the third example under Scope in the Interpolation documentation. Thanks for accepting my answer. By the way, a symbolic solution may be possible by using a double Integrate, instead of using ImplicitRegion. Jul 31, 2020 at 16:33

I now realize that the question can be solved analytically, for the most part, although not with ImplicitRegion. The constraints embodied in R1 and R2 can be solved to obtain θ in terms of v and parameters.

r1c1 = Reduce[θ (v + r) > r && v > 0, θ] // Last
(* θ > 1/(1 + 2 v) *)
r1c2 = Reduce[2 θ v s + (1 - θ) r > λ1 && v > 0 && s > 0 && s != 1/(4 v), θ] // Last
(* (0 < s < 1/(4 v) && θ < (-1 + 2 λ1)/(-1 + 4 s v)) ||
(s > 1/(4 v) && θ > (-1 + 2 λ1)/(-1 + 4 s v)) *)

r2c1 = Reduce[θ (v + r) < r && v > 0, θ] // Last
(* θ < 1/(1 + 2 v) *)
r2c2 = Reduce[2 θ v s + λ2 v > λ1 && v > 0 && s > 0, θ] // Last // Simplify
(* θ > (-0.196007 v + 0.5 λ1)/(s v) *)


The corresponding integrals then can be evaluated symbolically in just a few minutes, although the results are a bit long to reproduce here. (Edit: Apply N, if necessary, so that r2c2 has the form shown. Otherwise computations of f1i2 and f2i2 are very slow.)

f1i1 = Integrate[(1 - θ) r Boole[r1c1 && r1c2], {v, 0, 1}, {θ, 0, 1},
Assumptions -> λ1 > 0, GenerateConditions -> True];
f1i2 = Integrate[θ v Boole[r2c1 && r2c2], {v, 0, 1}, {θ, 0, 1},
Assumptions -> λ1 > 0 && s > 0, GenerateConditions -> True] // Simplify;
f2i1 = Integrate[Boole[r1c1 && r1c2], {v, 0, 1}, {θ, 0, 1},
Assumptions -> λ1 > 0, GenerateConditions -> True];
f2i2 = Integrate[Boole[r2c1 && r2c2], {v, 0, 1}, {θ, 0, 1},
Assumptions -> λ1 > 0 && s > 0, GenerateConditions -> True] // Simplify;


Determining s and λ1 for the illustrative conditions used in my earlier numerical solution yields

sols = FindRoot[{f1i1 + f1i2 == .03, f2i1 + f2i2 == .135}, {{s, .25}, {λ1, .35}}]
(* {s -> 0.144367, λ1 -> 0.356326} *)


differing from the results of the earlier entirely numerical calculations by about 1%.

A simpler code producing essentially the same results is, in its entirety,

r = 1/2;
λ2 = N[-(1/2) + r + 1/(1 + r) - r^2 Log[1 + 1/r]];

f1i1 = Integrate[(1 - θ) r Boole[θ (v + r) > r && 2 θ v s + (1 - θ) r > λ1],
{v, 0, 1}, {θ, 0, 1}, Assumptions -> λ1 > 0 && s > 0, GenerateConditions -> True];
f1i2 = Integrate[θ v Boole[θ (v + r) < r && 2 θ v s + λ2 v > λ1],
{v, 0, 1}, {θ, 0, 1}, Assumptions -> λ1 > 0 && s > 0, GenerateConditions -> True]
// Simplify;
f2i1 = Integrate[Boole[(1 - θ) r Boole[θ (v + r) > r && 2 θ v s + (1 - θ) r > λ1],
{v, 0, 1}, {θ, 0, 1}, Assumptions -> λ1 > 0 && s > 0, GenerateConditions -> True];
f2i2 = Integrate[Boole[θ (v + r) < r && 2 θ v s + λ2 v > λ1],
{v, 0, 1}, {θ, 0, 1}, Assumptions -> λ1 > 0 && s > 0, GenerateConditions -> True]
// Simplify;


after which s and λ1 can be determined as desired.