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} \begin{equation} \Omega = \{(x,y): (x/3-1)^2+y^2 \leq 1 \} \end{equation}
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.
