On Wed, 5 Apr 2006 11:10:14 +0000 (UTC), axman <yuerlong@gmail.com> wrote:
> Hi, Group
> I am new to Mathematica and now meet one problem I can't solve by
> self. First I defined 3 functions: h, f and g
>
> \!\(\(h[t_] := \(-20\)*Cos[t] Exp[\(-t\)\/20];\)\[IndentingNewLine]
> \(f[x_, y_] := Cos[x] Sin[y];\)\[IndentingNewLine]
> \(g[_, _, t_] :=
> If[t == 0, \(-1\), If[h[t] <
> , \(-1\), If[h[t] > , \(+1\), g[, , t - 0.1]]]];\)\)
>
>
> Now I try to integrate the function f[x,y]*g[x,y,10] on a triangle
> region:
>
> \!\(\_\(-20\)\%20\(\_\(-20\)\%x\((N[f[
> x, y]\ g[x, y, 10]])\) \[DifferentialD]y \[DifferentialD]x\)\)
>
> The output result is strange. How to get correct answer for this
> integration?
>
> Thank you for help.
>
This is a tricky one, because it involves recursion and so we can't just
unwind all the nested If's at once. One way to handle this is to introduce
a dummy function and perform the integrations step by step
(PiecewiseIntegrate code is available at
http://library.wolfram.com/infocenter/MathSource/5117/ ):
In[2]:= h[t_] := -20*Cos[t]*Exp[-t/20];
f[x_, y_] := Cos[x]*Sin[y];
g[a_, b_, t_] := If[t == 0, -1, If[h[t] < b, -1,
If[h[t] > a, 1, g[a, b, t - 1/10]]]];
In[5]:= g2[a_, b_, t_] := If[t == 0, -1, If[h[t] < b, -1,
If[h[t] > a, 1, g2dummy[a, b, t - 1/10]]]];
In[6]:= FixedPoint[# /.
{g2dummy -> g2, Integrate -> PiecewiseIntegrate}&,
Integrate[f[x, y]*g2[x, y, 10`20], {x, -20, 20}, {y, -20, x}]] //
Timing
Out[6]= {12.516*Second,
-11.4346710658068048816232977508`15.290914394480588}
To evaluate the integral by purely numerical methods we notice that the
problematic points are where x == h[t] or y == h[t], and t will take
values 10, 10 - 1/10, and so on:
In[7]:= Lcp = Union@ h@ Range[10., 0., -.1];
NIntegrate[f[x, y]*g[x, y, 10.],
{x, -20, Sequence @@ Lcp, 20} // Evaluate,
{y, -20, Sequence @@ (Min[#, x]&) /@ Lcp, x} // Evaluate,
Method -> MultiDimensional] // Timing
Out[8]= {48.406*Second, -11.434671066765237`}
In fact, with exact input we can even get the exact value of the integral
using the same combination of PiecewiseIntegrate and FixedPoint; the exact
answer is
(20*Cos[10])/Sqrt[E] - 2*Cos[20]*Sin[(20*Cos[31/10])/E^(31/200)]
+ 2*Cos[(20*Cos[31/5])/E^(31/100)]*Sin[(20*Cos[31/10])/E^(31/200)] -
2*Cos[(20*Cos[31/5])/E^(31/100)]*Sin[(20*Cos[47/5])/E^(47/100)]
+ 2*Cos[(20*Cos[10])/Sqrt[E]]*Sin[(20*Cos[47/5])/E^(47/100)] -
(1/2)*Sin[(40*Cos[10])/Sqrt[E]]
Maxim Rytin
m.r@inbox.ru
|