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: [mg80457] Re: Newbie question on FindRoot and NIntergrate
  • From: Jean-Marc Gulliet <jeanmarc.gulliet at gmail.com>
  • Date: Thu, 23 Aug 2007 01:00:19 -0400 (EDT)
  • Organization: The Open University, Milton Keynes, UK
  • References: <fagu3e$937$1@smc.vnet.net>

Hoa Bui wrote:

> 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] ?

FindRoot returns a list of transformation rules (transformation rules of 
the form y -> somenumber) which cannot be directly interpreted as a 
number by NIntegrate. Therefore, you must modified the definition of 
hh[x] to force it to return a number. To do so, we use the replacement 
operator /. in the definition below.

Also, hh[x] must be called only for numeric value of its argument x. 
This can be achieved by adding a constraint on x in the LHS of hh as 
done  below.

In[1]:= 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}};

In[2]:= linNx = Interpolation[points, InterpolationOrder -> 1];
ff[x_] := Integrate[linNx[s], {s, 0.0001, x}]

In[4]:= gg[x_] := 0.9 x^(1/0.9)

In[5]:= hh[x_?NumberQ] :=
   y /. FindRoot[ff[y] == gg[x], {y, 0.0001, 0.019}];

In[6]:= hh[0.001]

Out[6]= 0.00107654

In[7]:= hh[0.01]

Out[7]= 0.0101302

In[8]:= NIntegrate[hh[x], {x, 0.0001, 0.019}]

Out[8]= 0.000183073

-- 
Jean-Marc


  • Prev by Date: Re: Re: Redefine Arg to return a value from 0 to 2 pi
  • Next by Date: Re: RE: Re: Missing information in Mathematica 6.0 Documentation
  • Previous by thread: Re: Newbie question on FindRoot and NIntergrate
  • Next by thread: Re: Newbie question on FindRoot and NIntergrate