Re: Re: integration of piecewise convex bivariate function with 6 parameters
- To: mathgroup at smc.vnet.net
- Subject: [mg79074] Re: [mg78965] Re: integration of piecewise convex bivariate function with 6 parameters
- From: "Michael Chen" <vancouver.michael at gmail.com>
- Date: Wed, 18 Jul 2007 02:50:37 -0400 (EDT)
- References: <f74t0t$91f$1@smc.vnet.net> <200707131018.GAA12804@smc.vnet.net>
Thanks Maxim, I cannot follow the installation instruction on my version 6
though. There is no "File->Save as Special->Package Format". I did copied
the code as pure text, but have no clue what to do next. Would you please
give more details? Thanks.
Michael Chen
On 7/13/07, m.r at inbox.ru <m.r at inbox.ru> wrote:
>
> On Jul 12, 4:41 am, "Michael Chen" <vancouver.mich... at gmail.com>
> wrote:
> > Dear there,
> >
> > I spent roughly one week to manually integrate a piecewise convex
> function
> > with 6 parameters, and I can't imagine how happy I would be if
> Mathematica
> > can do it for me too. Following is the function, and the Mathematica
> command
> > I tried.
> >
> > f(h1,h1)= if (h1 - x1 <= 0 or h2 - x2 <= 0) then | h1 - x1| + | h2 - x2
> |
> > else max( h1 - x1, h2 - x2), where x1, x2 are real parameters. And I
> > integrate the function f(h1, h2) over rectangular region {h1, B1, U1},
> {h2,
> > B2, U2}, where B1, U1, B2, U2 are reals, and B1 < U1, B2 < U2.
> >
> > I tried this in Mathematica 6.0:
> >
> > g[x1_Real, x2_Real] := Integrate[
> > If [h1 - x1 <= 0 || h2 - x2 <= 0, abs[h1 - x1] + abs[h2 - x2] ,
> > max[h1 - x1, h2 - x2]], {h1, B1, U1}, {h2, B2, U2},
> > Assumptions ->
> > Im[B1] == 0 && Im[B2] == 0 && Im[U1] == 0 && Im[U2] == 0 &&
> > Im[h1] == 0 && Im[h2] == 0 && B1 < U1 && B2 < U2 ];
> >
> > However there is simply no response from Mathematica 6.0 after I press
> > Shift+Enter. Could anybody give me any suggestions? Thanks a lot.
> >
> > ( The example
> onhttp://reference.wolfram.com/mathematica/tutorial/IntegralsOverRegion...,
> > [image: Click for copyable input] <javascript:input('i_9')>, seems very
> > close to what I am doing. And Mathematica gives back a piecewise
> function
> > in a as answer. I know the integral of my function f(h1,h2) is a
> > piecewise function of (x1,x2) too.)
> >
> > Michael Chen
>
> Things that generally help:
>
> capitalize all the function names (abs -> Abs, max -> Max);
>
> make the change of variables to reduce the number of parameters (x1
> and x2 can be substituted later);
>
> perform the integration step by step (integrating out one variable at
> a time) and term by term;
>
> use PiecewiseIntegrate (
> http://library.wolfram.com/infocenter/MathSource/5117/
> ).
>
> In[2]:= f[h1_, h2_] := If[h1 <= 0 || h2 <= 0,
> Abs[h1] + Abs[h2], Max[h1, h2]]
>
> In[3]:= ii = PiecewiseIntegrate[f[h1, h2], {h1, B1, U1},
> Assumptions -> B1 < U1]
>
> Out[3]= If[0 < B1 && h2 <= 0, 1/2 (-B1^2 + 2 B1 h2 - 2 h2 U1 + U1^2),
> 0] +
> If[B1 < 0 && U1 < 0, 1/2 (B1^2 - U1^2 - 2 B1 Abs[h2] + 2 U1 Abs[h2]),
> 0] + If[B1 < 0 && 0 <= U1, 1/2 (B1^2 - 2 B1 Abs[h2]), 0] +
> If[B1 == 0 && h2 <= 0 && 0 < U1, 1/2 (-2 h2 U1 + U1^2), 0] +
> If[0 < B1 && h2 == B1 && h2 < U1, 1/2 (-h2^2 + U1^2), 0] +
> If[0 < B1 && 0 < h2 && h2 < B1, 1/2 (-B1^2 + U1^2), 0] +
> If[0 < B1 && B1 < h2 && U1 == h2, -h2 (B1 - U1), 0] +
> If[0 < B1 && B1 < h2 && h2 < U1, 1/2 (-2 B1 h2 + h2^2 + U1^2), 0] +
> If[0 < B1 && B1 < h2 && U1 < h2, -h2 (B1 - U1), 0] +
> If[B1 < 0 && h2 <= 0 && 0 < U1, 1/2 (-2 h2 U1 + U1^2), 0] +
> If[B1 <= 0 && 0 < h2 && U1 == h2, h2 U1, 0] +
> If[B1 <= 0 && 0 < h2 && h2 < U1, 1/2 (h2^2 + U1^2), 0] +
> If[B1 <= 0 && 0 < h2 && 0 < U1 && U1 < h2, h2 U1, 0]
>
> In[4]:= ii = PiecewiseIntegrate[#, {h2, B2, U2},
> Assumptions -> B2 < U2]& /@ ii
>
> Out[4]= If[0 < B1 && B1 < B2 && U1 < B2,
> 1/2 (B1 B2^2 - B2^2 U1 - B1 U2^2 + U1 U2^2), 0] +
> If[0 < B1 && B2 < 0 && 0 < U2,
> 1/2 (B1^2 B2 - B1 B2^2 + B2^2 U1 - B2 U1^2), 0] +
> If[0 < B1 && B2 < 0 && U2 <= 0,
> 1/2 (B1^2 B2 - B1 B2^2 + B2^2 U1 - B2 U1^2 - B1^2 U2 + U1^2 U2 +
> B1 U2^2 - U1 U2^2), 0] +
> If[0 < B1 && B2 <= 0 && U2 == B1, -1/2 (B1^2 - U1^2) U2, 0] +
> If[0 < B1 && B2 <= 0 && B1 < U2, -1/2 B1 (B1^2 - U1^2), 0] +
> If[B1 < 0 && 0 <= B2 && U1 < 0,
> 1/2 (-B1^2 B2 + B1 B2^2 - B2^2 U1 + B2 U1^2 + B1^2 U2 - U1^2 U2 -
> B1 U2^2 + U1 U2^2), 0] +
> If[B1 < 0 && 0 <= B2 && 0 <= U1,
> 1/2 (-B1^2 B2 + B1 B2^2 + B1^2 U2 - B1 U2^2), 0] +
> If[B1 == 0 && B2 < 0 && 0 < U1 && 0 < U2, 1/2 (B2^2 U1 - B2 U1^2),
> 0] + If[B1 == 0 && B2 < 0 && 0 < U1 && U2 <= 0,
> 1/2 (B2^2 U1 - B2 U1^2 + U1^2 U2 - U1 U2^2), 0] +
> If[0 < B1 && B2 == B1 && U1 < B1 && B1 < U2,
> 1/2 (B2^3 - B2^2 U1 - B2 U2^2 + U1 U2^2), 0] +
> If[0 < B1 && B2 == B1 && B1 <= U1 && U1 < U2,
> 1/2 (B2 U1^2 - U1^3 - B2 U2^2 + U1 U2^2), 0] +
> If[0 < B1 && 0 < B2 && B2 < B1 && U2 == B1,
> 1/2 (B1^2 - U1^2) (B2 - U2), 0] +
> If[0 < B1 && 0 < B2 && B2 < B1 &&
> B1 < U2, -1/2 (B1 - B2) (B1^2 - U1^2), 0] +
> If[0 < B1 && 0 < B2 && B2 < B1 && U2 < B1,
> 1/2 (B1^2 - U1^2) (B2 - U2), 0] +
> If[0 < B1 && B1 < B2 && B2 < U1 && U2 == U1,
> 1/6 (3 B1 B2^2 - B2^3 - 3 B1 U2^2 - 3 B2 U2^2 + 4 U2^3), 0] +
> If[0 < B1 && B1 < B2 && B2 < U1 && U1 < U2,
> 1/6 (3 B1 B2^2 - B2^3 - 3 B1 U1^2 - 3 B2 U1^2 + 4 U1^3), 0] +
> If[0 < B1 && B1 < B2 && B2 < U1 && U2 < U1,
> 1/6 (3 B1 B2^2 - B2^3 - 3 B2 U1^2 + 3 U1^2 U2 - 3 B1 U2^2 + U2^3),
> 0] + If[0 < B1 && B1 < B2 && B2 <= U1 && U1 < U2,
> 1/2 (B1 U1^2 - U1^3 - B1 U2^2 + U1 U2^2), 0] +
> If[0 < B1 && B2 < B1 && B1 < U1 && U1 < U2,
> 1/2 (B1 U1^2 - U1^3 - B1 U2^2 + U1 U2^2), 0] +
> If[0 < B1 && B2 < B1 && U1 <= B1 && B1 < U2,
> 1/2 (B1^3 - B1^2 U1 - B1 U2^2 + U1 U2^2), 0] +
> If[0 < B1 && B2 <= 0 && 0 < U2 && U2 < B1, -1/2 (B1^2 - U1^2) U2,
> 0] + If[0 < B1 && B2 <= B1 && B1 < U1 && U2 == U1,
> 1/3 (B1^3 - 3 B1 U2^2 + 2 U2^3), 0] +
> If[0 < B1 && B2 <= B1 && B1 < U1 && U1 < U2,
> 1/3 (B1^3 - 3 B1 U1^2 + 2 U1^3), 0] +
> If[B1 < 0 && B2 < 0 && 0 < U1 && 0 < U2, 1/2 (B2^2 U1 - B2 U1^2),
> 0] + If[B1 < 0 && B2 < 0 && 0 < U1 && U2 <= 0,
> 1/2 (B2^2 U1 - B2 U1^2 + U1^2 U2 - U1 U2^2), 0] +
> If[B1 < 0 && B2 < 0 && U1 < 0 && 0 < U2,
> 1/2 (B1^2 U2 - U1^2 U2 - B1 U2^2 + U1 U2^2), 0] +
> If[B1 < 0 && B2 < 0 && U1 < 0 && U2 < 0,
> 1/2 (-B1^2 B2 - B1 B2^2 + B2^2 U1 + B2 U1^2 + B1^2 U2 - U1^2 U2 +
> B1 U2^2 - U1 U2^2), 0] +
> If[B1 < 0 && B2 < 0 && U1 < 0 && 0 <= U2,
> 1/2 (-B1^2 B2 - B1 B2^2 + B2^2 U1 + B2 U1^2), 0] +
> If[B1 < 0 && B2 < 0 && 0 <= U1 && 0 < U2, 1/2 (B1^2 U2 - B1 U2^2),
> 0] + If[B1 < 0 && B2 < 0 && 0 <= U1 && U2 < 0,
> 1/2 (-B1^2 B2 - B1 B2^2 + B1^2 U2 + B1 U2^2), 0] +
> If[B1 < 0 && B2 < 0 && 0 <= U1 && 0 <= U2, 1/2 (-B1^2 B2 - B1 B2^2),
> 0] + If[B1 <= 0 && 0 < B2 && 0 < U1 && U1 < B2,
> 1/2 (-B2^2 U1 + U1 U2^2), 0] +
> If[B1 <= 0 && 0 < B2 && B2 < U1 && U2 == U1,
> 1/6 (-B2^3 - 3 B2 U2^2 + 4 U2^3), 0] +
> If[B1 <= 0 && 0 < B2 && B2 < U1 && U1 < U2,
> 1/6 (-B2^3 - 3 B2 U1^2 + 4 U1^3), 0] +
> If[B1 <= 0 && 0 < B2 && B2 < U1 && U2 < U1,
> 1/6 (-B2^3 - 3 B2 U1^2 + 3 U1^2 U2 + U2^3), 0] +
> If[B1 <= 0 && 0 < B2 && B2 <= U1 && U1 < U2, 1/2 (-U1^3 + U1 U2^2),
> 0] + If[B1 <= 0 && B2 <= 0 && 0 < U1 && U2 == U1, (2 U2^3)/3, 0] +
> If[B1 <= 0 && B2 <= 0 && 0 < U1 && U1 < U2, (2 U1^3)/3, 0] +
> If[B1 <= 0 && B2 <= 0 && 0 < U1 && U1 < U2, 1/2 (-U1^3 + U1 U2^2),
> 0] + If[0 < B1 && B2 <= B1 && B1 < U1 && B1 < U2 && U2 < U1,
> 1/6 (2 B1^3 - 3 B1 U1^2 + 3 U1^2 U2 - 3 B1 U2^2 + U2^3), 0] +
> If[B1 <= 0 && B2 <= 0 && 0 < U1 && 0 < U2 && U2 < U1,
> 1/6 (3 U1^2 U2 + U2^3), 0]
>
> In[5]:= ii = ii /. {p1: B1 | U1 :> p1 - x1, p2: B2 | U2 :> p2 - x2};
>
> A numeric check:
>
> In[6]:= Table[(ii - NIntegrate[f[h1 - x1, h2 - x2],
> {h1, B1, U1}, {h2, B2, U2}]) Boole[B1 < U1 && B2 < U2] /.
> Thread[{B1, U1, B2, U2, x1, x2} -> RandomReal[{-10, 10}, 6]],
> {100}] // Abs // Max
>
> Out[6]= 9.094947*10^-13
>
> Maxim Rytin
> m.r at inbox.ru
>
>
>
- References:
- Re: integration of piecewise convex bivariate function with 6 parameters
- From: m.r@inbox.ru
- Re: integration of piecewise convex bivariate function with 6 parameters