1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Mauricio, How about doing a Laplace Transformation of the equations and solve in the Laplace Domain? $\endgroup$ – Zviovich Oct 1 '15 at 14:31
  • $\begingroup$ @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. $\endgroup$ – Mauricio Fernández Oct 1 '15 at 14:42
3
$\begingroup$
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[] :)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.