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.
