MathGroup Archive 2007

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

Search the Archive

Re: Newbie question on FindRoot and NIntergrate

  • To: mathgroup at smc.vnet.net
  • Subject: [mg80489] Re: [mg80437] Newbie question on FindRoot and NIntergrate
  • From: DrMajorBob <drmajorbob at bigfoot.com>
  • Date: Thu, 23 Aug 2007 01:16:54 -0400 (EDT)
  • References: <22117982.1187821637705.JavaMail.root@m35>
  • Reply-to: drmajorbob at bigfoot.com

Your problem is quite a tangle (IMHO), so forgive me if I've gotten  
offtrack somewhere.

I start where you did:

points = {{0.0001, 0.359381}, {0.0002866, 0.403984}, {0.000821394,
     0.454122}, {0.00235411, 0.510482}, {0.00674688,
     0.573838}, {0.0193365, 0.645056}};
linNx = Interpolation[points, InterpolationOrder -> 1];
{xmin, xmax} = points[[{1, -1}, 1]];
Plot[linNx[x], {x, xmin, xmax}]

Your ff I call f, and I compute it using NDSolve, not Integrate:

Clear[f, x]
f[x_] = f[x] /.
    First@NDSolve[{f'[x] == linNx[x], f[xmin] == 0},
      f[x], {x, xmin, xmax}];
Plot[f[x], {x, xmin, xmax}]

If fi is the inverse function to f, then y so that f[y] == g[x] is just y  
= fi@g@x.

To find fi, we'll need the upper and lower values of f, and we need to  
know that f is monotone increasing:

{fmin, fmax} = f /@ {xmin, xmax}

{4.81833*10^-23, 0.0110942}

{fmin, fmax} = f /@ {xmin, xmax}
f'[xmin]
Plot[f'[x], {x, xmin, xmax}]

{4.81833*10^-23, 0.0110942}

0.359381

f'[x] > 0 across the domain, so we can find the inverse using NDSolve  
again:

Clear[fi, y]
fi[y_] = fi[y] /.
    First@NDSolve[{linNx[fi[y]] fi'[y] == 1, fi[fmin] == xmin},
      fi, {y, fmin, fmax}];
Plot[fi[y], {y, fmin, fmax}]

Now define g and h (your gg and hh):

Clear[g, x, h]
g[x_] = 0.9 x^(1/0.9);
h[x_] = fi@g@x;

and find the upper and lower limits for the domain of h. That's where g 
takes on the values fmax and fmin:

Off[Solve::"ifun"]
{gmin, gmax} = Flatten[x /. Solve[g[x] == #, x] & /@ {fmin, fmax}]

{9.03209*10^-21, 0.0191323}

h[x] is almost exactly x,:

Plot[h[x] - x, {x, gmin, gmax}, PlotRange -> All]

(very small values) and that's because f and g are almost exactly the same:

Plot[f@x - g@x, {x, xmin, xmax}]

so fi is very close to an inverse of g.

Finally, integrate h over its domain:

NIntegrate[h[x], {x, gmin, gmax}]

0.000185636

Bobby

On Wed, 22 Aug 2007 03:49:45 -0500, Hoa Bui <hoabui05 at gmail.com> wrote:

> Dear all,
>
> Please help me with this problem:
>
> I have a series of points:
> In[9]:=points
> Out[9]={{0.0001,0.359381},{0.0002866,0.403984},{0.000821394,0.454122},{0.00235411,0.510482},{0.00674688,0.573838},{0.0193365,0.645056}}
>
> Define one of my function as the trapezoidal area created by these
> points up to some x:
> linNx = Interpolation[points, InterpolationOrder -> 1];
> ff[x_] := Integrate[linNx[s], {s, 0.0001, x}]
>
> Define the second function as a power law:
> gg[x_] := 0.9 x^(1/0.9)
>
> My third function is the root of the equation:
> hh[x_] := FindRoot[ff[y] == gg[x], {y, 0.0001, 0.019}];
>
> I can evaluate h at specific values of x, e.g. h[0.001] ({y ->
> 0.00107654}), h[0.01] ({y -> 0.0101302}), etc...
>
> Now say I want to use hh[x] in an integral,
> NIntergrate[hh[x],{x, 0.0001, 0.019}]
> obviously it doesn't work, and I have tried using Solve instead of
> FindRoot but it also did not output a numerical value for the integral
> because the inverse function is not in closed form or something..
>
> Is there a way for me to compute the integral of hh[x] ?
>
> Thank you all so much,
> Hoa Bui
>
>



-- 

DrMajorBob at bigfoot.com


  • Prev by Date: Re: Re: Intensive numerical calculations
  • Next by Date: Re: Re: Re: Cell Bracket Symbols
  • Previous by thread: Re: Newbie question on FindRoot and NIntergrate
  • Next by thread: ListSurfacePlot3D in Mathematica Version 6