How to solve nonlinear equations involving integrals - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-10-15T22:11:54Z https://mathematica.stackexchange.com/feeds/question/95931 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/95931 1 How to solve nonlinear equations involving integrals Mauricio Fernández https://mathematica.stackexchange.com/users/24710 2015-10-01T14:14:48Z 2015-10-09T14:45:13Z <p>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: <a href="https://mathematica.stackexchange.com/questions/91333/how-to-solve-algebra-equations-containing-integration-and-parameters">How to solve algebra equations containing containing integration and parameters</a> . 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 <strong>constant</strong> 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.</p> <p>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 &amp;=&amp; 16.9381 \\ \int_\Omega \exp(a_1 \sin(x)y)\cos(a_2 x) dx dy &amp;=&amp; -21.057 \end{eqnarray} \begin{equation} \Omega = \{(x,y): (x/3-1)^2+y^2 \leq 1 \} \end{equation}</p> <p>QUESTION: ANY IDEAS ON HOW TO TREAT SUCH A PROBLEM WITH MATHEMATICA FUNCTIONS? </p> <p>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. </p> <p>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.</p> <p><strong>Minimal example in 2D</strong></p> <p>Integration region, parametrized integral and right hand side</p> <pre><code>(*Integration region*) reg = ImplicitRegion[(x/3 - 1)^2 + (y + 0)^2 &lt;= 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 </code></pre> <p>Solution with FindRoot and Integrate</p> <pre><code>(*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 </code></pre> <p>Solution with rudimentary Newton algorithm</p> <pre><code>(*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 -&gt; a1, a2s -&gt; 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] &gt; tol &amp;&amp; counter &lt; 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 </code></pre> <p>Thank you very much for any idea or information about a duplicate (sorry in that case).</p> <p><strong>edit 2015-Oct-02:</strong> 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.</p> https://mathematica.stackexchange.com/questions/95931/-/96017#96017 3 Answer by Mariusz Iwaniuk for How to solve nonlinear equations involving integrals Mariusz Iwaniuk https://mathematica.stackexchange.com/users/26828 2015-10-02T11:24:49Z 2015-10-02T11:34:45Z <pre><code>Clear["Global`*"] reg = ImplicitRegion[(x/3 - 1)^2 + (y + 0)^2 &lt;= 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 -&gt; 20] </code></pre> <blockquote> <p>$\{\text{a1}\to 7.0000001153135245011,\text{a2}\to 1.0471978494085482801\}$</p> </blockquote> <p><strong>Check solution:</strong></p> <pre><code>{int1[a1 /. sol, a2 /. sol], int2[a1 /. sol, a2 /. sol]} </code></pre> <blockquote> <p>$\{16.9381,-21.057\}$</p> </blockquote> <p>You can change, a starting points <code>(a1 ,a2)</code> in <strong>FindRoot[] :)</strong></p>