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