How to solve nonlinear equations involving integrals

Sorry, if this is a duplicate, I was not able to find a corresponding answer. The only question related to mine, as far as I can see it, is this one: How to solve algebra equations containing containing integration and parameters . But no answers were given. My problem: I have a nonlinear algebraic system for some unknown variables, say, $a_1$ and $a_2$, which are in a nonlinear integral of a function containing $a_1$ and $a_2$ over a region $\Omega$, e.g., an ellipse. The constant right hand side is also given. I really just have a nonlinear system. The integrals are not elementary, so I have to use some trick.

Example problem in 2D, unknowns are $a_1$ and $a_2$ \begin{eqnarray} \int_\Omega \exp(a_1 \sin(x)y)\sin(a_2 x) dx dy &=& 16.9381 \\ \int_\Omega \exp(a_1 \sin(x)y)\cos(a_2 x) dx dy &=& -21.057 \end{eqnarray} $$\Omega = \{(x,y): (x/3-1)^2+y^2 \leq 1 \}$$

QUESTION: ANY IDEAS ON HOW TO TREAT SUCH A PROBLEM WITH MATHEMATICA FUNCTIONS?

Until now, I am able to treat this using FindRoot and Integrate, but already for the 2D example given below this is kind of slow. I dont know if there is a way to use NIntegrate there. Alternatively, I tried a rudimentary Newton algorithm and seems ok, see code below.

I am aware that I could use some quadrature rule in order to "overcome" the integrals and formulate a purely algebraic system. Sadly, later I have to treat a high dimensional problem over a complicated set and the quadrature rules are not very helpful. I am also aware that I could use Monte Carlo integration. I tried this but in the set I have to solve my problem later the results just dont get better and at some point I just run out of memory.

Minimal example in 2D

Integration region, parametrized integral and right hand side

(*Integration region*)
reg = ImplicitRegion[(x/3 - 1)^2 + (y + 0)^2 <= 1, {x, y}];
(*Numerical evaluation of parametrized integral*)
integrand[x_, y_, a1_, a2_] := {
1.*Exp[a1*Sin[x]*y]*Sin[a2*x],
1.*Exp[a1*Sin[x]*y]*Cos[a2*x]
};
Nint[a1_, a2_] :=NIntegrate[integrand[x, y, a1, a2], Element[{x, y}, reg]] // Quiet;
(*Right hand side*)
atest = {7, Pi/3} // N
rhs = Nint@@atest


Solution with FindRoot and Integrate

(*Solve with Mathematica function FindRoot and symbolic integration*)
start = DateString[]
root = FindRoot[
Integrate[integrand[x, y, a1, a2], Element[{x, y}, reg]] -
rhs, {{a1, 5}, {a2, Pi/4}}]
end = DateString[]
DateDifference[start, end, {"Hour", "Minute", "Second"}]
aroot = {a1, a2} /. root;
Nint@@aroot


Solution with rudimentary Newton algorithm

(*Rudimentary Newton algorithm, based on Mathematica functions*)
(*Numerical evaluation of Jacobian of parametrized integral*)
NintJac[a1_, a2_] :=
NIntegrate[
D[integrand[x, y, a1s, a2s], {{a1s, a2s}, 1}] /. {a1s -> a1,
a2s -> a2}, Element[{x, y}, reg]] // Quiet;
(*Newton settings*)
alast = {5, Pi/4};(*first guess*)
tol = 10^(-5);
counter = 0;
error = rhs - Nint@@alast;
maxit = 10^2;
information = {counter, alast // N, Norm[error] // ScientificForm};
Print["Start with"];
Print[information];
(*Newton*)
start = DateString[]
Monitor[
While[
Norm[error] > tol && counter < maxit
,
counter = counter + 1;
Jloc = NintJac@@alast;
linsol = LinearSolve[Jloc, error];
alast = alast + linsol;
error = rhs - Nint@@alast;
information = {counter, alast // N, Norm[error] // ScientificForm};
]
, information
]
end = DateString[]
DateDifference[start, end, {"Hour", "Minute", "Second"}]
(*Print results*)
information


Thank you very much for any idea or information about a duplicate (sorry in that case).

edit 2015-Oct-02: sorry, I should have given the equations from the beginning in Latex. See now example problem in 2D at the beginning. It's the same example as in the code.

• Mauricio, How about doing a Laplace Transformation of the equations and solve in the Laplace Domain? – Zviovich Oct 1 '15 at 14:31
• @Zviovich hmmm... can you explain a little? I know that the laplace transformation can be applied to differential equations in order to solve them algebraically and then transform them back. But how does that apply to this problem? Sorry, I dont see it. – Mauricio Fernández Oct 1 '15 at 14:42

1 Answer

Clear["Global*"]

reg = ImplicitRegion[(x/3 - 1)^2 + (y + 0)^2 <= 1, {x, y}];

int1[a1_?NumericQ, a2_?NumericQ] :=
NIntegrate[Exp[a1*Sin[x]*y]*Sin[a2*x], {x, y} \[Element] reg];
int2[a1_?NumericQ, a2_?NumericQ] :=
NIntegrate[Exp[a1*Sin[x]*y]*Cos[a2*x], {x, y} \[Element] reg];

sol = FindRoot[{int1[a1, a2] == Rationalize[16.9381, 0],
int2[a1, a2] == Rationalize[-21.057, 0]}, {{a1, 7}, {a2, Pi/3}},
WorkingPrecision -> 20]


$\{\text{a1}\to 7.0000001153135245011,\text{a2}\to 1.0471978494085482801\}$

Check solution:

{int1[a1 /. sol, a2 /. sol], int2[a1 /. sol, a2 /. sol]}


$\{16.9381,-21.057\}$

You can change, a starting points (a1 ,a2)` in FindRoot[] :)