Re: definite numerical integration

[David Bateman]

roberto wrote:

> On 9/21/05, David Bateman <address@hidden> wrote:
>  
>
> > Then all you can do is something like
> >
> > x0 = ??;
> > y0 =  f(x0);
> > a = ??;
> > b = ??;
> >
> > # x0 must be montonically increasing
> > [x0, idx] = sort(x0);
> > y(idx) = y0;
> >
> > # Treat core of the integral with simpson's rule
> > Istart = find (x0 >= a)(1);
> > Istop = find (x0 <= b)(end);
> > Int = sum(0.5 * (y0((Istart+1):Istop) + y0(Istart:(Istop-1))) .*
> > (x0((Istart+1):Istop) - x0(Istart:(Istop-1))));
> >
> > # Treat the tail
> > if (x0(Istart) != a)
> >  Int += 0.5 * (x0(Istart) - a) * (y0(Istart) + y0(Istart-1));
> > endif
> > if (x0(Istop) != b)
> >  Int += 0.5 * (b - x0(Istop)) * (y0(Istop+1) + y0(Istop));
> > endif
> > D.
> >    
> well thank you for contribution then;
> i'll ask one more point: my function to be integrated has also a parameter "c";
> now i want to solve the integral and then find out a value for this
> parameter such that the value of the integral when the upper bound of
> integration is e.g. "r_s" should be equal to a given value "t_s";
> is it possible to solve definite integral of functions depending also
> on a parameter to be fixed to satisfy a condition?
>
> thank you very much again,
> if the problem is known, please address me to some reference
>
>  
I'm not really sure I understand what you won't, but think its a formula 
with a parameter c which is unknown that you want the definite integral 
over some range to give a certain value and from that solve for c. I 
think you should be able to you something like the above with the 
addition of an optimization loop. Something like the below will probably 
do it, though I haven't bother testing and the only octave parser its 
gone through is the one in my brain..

D.

function Int = myquad (x0, y0, a, b)

 # x0 must be montonically increasing
 [x0, idx] = sort(x0);
 y(idx) = y0;

 # Treat core of the integral with simpson's rule
 Istart = find (x0 >= a)(1);
 Istop = find (x0 <= b)(end);
 Int = sum(0.5 * (y0((Istart+1):Istop) + y0(Istart:(Istop-1))) .*
 (x0((Istart+1):Istop) - x0(Istart:(Istop-1))));

 # Treat the tail
 if (x0(Istart) != a)
    Int += 0.5 * (x0(Istart) - a) * (y0(Istart) + y0(Istart-1));
 endif
 if (x0(Istop) != b)
   Int += 0.5 * (b - x0(Istop)) * (y0(Istop+1) + y0(Istop));
 endif
endfunction

function v = fun (c)
 global a b x0 v0;
 y0 = f(c, x0);  ## This is your function
 v = myquad (x0, y0, a, b) - v0;
endfunction

global a b x0 v0;

x0 = ...;   ## Your abssicca
c0 = ...;   ## Initial value of c
v0 = ...;   ## Desired value of the integral
a =  ...;   ## lower limit of definite integral
b = ...;   ## upper limit of definite integral

c = fsolve ('fun', c0);


--
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------