# 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? – bbgodfrey Jul 31 '20 at 3:12
• Also, supposing that Integrate can be performed, how are s and λ1 to be determined from the result? – bbgodfrey Jul 31 '20 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. – user391830 Jul 31 '20 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. – Ulrich Neumann Jul 31 '20 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. – bbgodfrey Jul 31 '20 at 13:47
• Thanks it' s very helpful! – Ulrich Neumann Jul 31 '20 at 13:59
• Thank you, this is super helpful! I have a syntax question, why do you need to Flatten the table? – user391830 Jul 31 '20 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. – bbgodfrey Jul 31 '20 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.