tom

Hi,
I need to integrate a function f(x) twice numerically, and while I
first though it should be rather straightforward, I am now doubting;
note that I am *not* talking of multidimensional integration, but
rather of something of this form:
G(x)=\int [ \int f(x) dx ] dx
where unfortunately it is not possible to provide the antiderivative of
f(x) in closed form. For this reason, it seems to me that it is not
possible to apply Gauss-type (or any other of the usual) methods to the
integral by a recursive call to a function, as I had thought in the
first instance.
The only way I can think of is to express f(x) as a power series and
integrate that twice up to some term, but I wonder whether this is a
accurate and a good way of doing it.
Surely there must be better methods - any ideas? Or have I
misunderstood something more basic in the numerical integration method
I thought of initially?
Tom

spellucci

In article <1142008361.919151.34120@p10g2000cwp.googlegroups. com>,
"Tom" <flurboglarf@mailinator.com> writes:


no, you can do it numerically straightforward, but you are not quite clear in
your notation: do you mean a table of the second antiderivative or a definite
integral?
some hints:


G(x)= \int_{a to x} \int_{a to z} f(t)dt dz = \int_{a to x} (x-t)f(t)dt
^^^ F(z) ^^^^^^^^^^^^
F(z) being the first antiderivative of f.

for other, similar cases, you could do it recursively: first tabulate
the first antiderivative using adaptive quadrature with a sufficiently
fine table of fixed nodes. then use this , combined with interpolation
in that table, for an adaptive quadrature in the outer integral.
hth
peter

dave dodson

How about using an ODE solver to integrate the system

y1' = f(x)
y2' = y1

Dave

One way: use a numerical integrator for the ODE y"=f(x).
For example, half station central differences (also known by many other
names). Another way: integrate the moment equation.

tom

Thanks to Peter (and also to the others who replied). That hint was
what I needed for this first step. However, that was only a preliminary
try for something uglier, namely things like this:
H(x)=\int_{a to x} [ g(z) \int_{a to z} f(t) dt ] dz
or
H(x)=\int_{a to x} exp[ \int_{a to z} f(t) dt ] dz
I had hoped that a solution for such cases would follow from the easier
problem, but it doesn't seem so (or does it?). Whereas it seems that
for the first equation one might be able to do a trick with partial
integration, I don't see possibilities for the second.
Is there any hope for a nice straightforward thing as well, or is the
combination of tabulation and interpolation the only way to do it?
I would rather avoid ODE solvers for this.
Tom

dave dodson

Tom writes:

Why do you want to avoid using an ODE solver? It seems like the easiest
and most straightforward way. For the first problem above, you can
solve the system:

y1' = f(x)
y2' = g(x) * y1

and for the second system, you can solve:

y1' = f(x)
y2' = exp(y1)

Dave

fred chapman

Tom, am I correct in assuming that you want an easily computed
*formula* which will approximate your function H(x)? Is that why you
want to avoid ODE solvers -- because they provide a table of values
rather than a formula? If the answer to both questions is yes, then
your problem is really a problem in approximation theory.

Some details would help us to help you better:

1. For what values of x do you want to approximate H(x)?

2. How accurately do you want to approximate H(x)?

3. What are typical functions g and f in your formulas for H(x)?

4. What do you want to do with the approximation for H(x)?

--
Frederick W. Chapman, Postdoctoral Fellow, University of Waterloo
http://www.scg.uwaterloo.ca/~fwchapman/

If you want to avoid ODE solvers, use the moment method: the
general solution of y"=f(x) with ICs y(0)=y0, y'(0)=y0p is

y(x) = y0 + y0p*x + int_0^x (x-y) f(y) dy

which involves only one integral, the convolution x*f, which
can be done numerically for each x.

fred krogh

If you don't want to work too hard at this, I believe the
multidimensional software at
http://mathalacarte.com/cb/mom.fcg/ya63 may be what you are looking for.
It is I believe the most efficient software out there for
multidimensional integrals up to dimension 3 or 4 or so, and will let
you take advantage of some of the structure that is in your problem.

If you want to work a bit harder, you may save computing time by making
up a new function, say q(z) = \int_a^z f(t) dt. This could be done by
doing an indefinite integration. Ideally there would be quadrature
software for doing this, but I don't know of it. The ODE software at
http://mathalacarte.com/cb/mom.fcg/ya64 has provision for saving the
solution, and thus you could save the solution to y' = f(t) over the
interval [a, x] using this software and then use interpolation (provided
by this software) as the definition of q(z). The former approach would
require considerably less work on your part. If the second approach is
used, you can save a bit of work by skipping the evaluation of f(t)
after the corrected value is obtained.

If you actually want H(x) for a bunch of different x values which can be
obtained in sequential order, you may find it faster to use an ode
solver to get these once you have a q(z).

Unless you are doing a lot of these computations with functions f and g
that are expensive to compute the quadrature software is what I would
recommend. Also if f(t) is much more expensive to compute that g(z) (in
your first case) you would probably be ahead to use integration by parts
to make the inner integral be an integral of g.

Note, I am the vendor of this software and depending on your use, you
may be expected to pay for it if you use it. Details are on the site at
http://mathalacarte.com/.

Regards,
Fred

tom

Thanks for the lot of suggestions, especially from Dave Dodson and Fred
Chapman. As it is apparently not possible to do it with a simple, neat
single integral, I'll continue entertaining myself trying out both the
interpolation+tabulation and the ODE thing.
As for Fred's questions: the purpose is really to express integrals in
the general solution of an ODE, which cannot be represented in closed
form for some special cases I'm looking at. To have as precise
solutions as possible, I had intended to calculate those terms of the
solution which can be done in closed form by programming the stuff
directly and do the integrals (which contain ugly combinations of exp
and erf) numerically.
Thanks again,
Tom