MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Integrate problem

  • To: mathgroup at smc.vnet.net
  • Subject: [mg65541] Re: Integrate problem
  • From: Maxim <m.r at inbox.ru>
  • Date: Thu, 6 Apr 2006 06:53:20 -0400 (EDT)
  • References: <e108im$lgb$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

On Wed, 5 Apr 2006 11:10:14 +0000 (UTC), axman <yuerlong at 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 at inbox.ru


  • Prev by Date: Re: graphic inside a gridbox
  • Next by Date: Re: ListDensityPlot and GraphicsArray
  • Previous by thread: Integrate problem
  • Next by thread: another problem with Infinite Products