[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: definite numerical integration
From: |
David Bateman |
Subject: |
Re: definite numerical integration |
Date: |
Thu, 22 Sep 2005 09:53:36 +0200 |
User-agent: |
Mozilla Thunderbird 0.8 (X11/20040923) |
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
-------------------------------------------------------------