Matlab - Unevenly spaced discretization for integration

Author

Unevenly spaced discretization for integration

neurino

2005-12-09, 7:26 pm

Hallo everybody.
I'm doomed by doubts.

Let's assume i need to solve an integral of product of a vector K.
The integrand oscillates dramatically at far right end.
I'm *not* interested in these high frequencies, then i sample K with an
unevenly spaced points
(that is more far to the right, bigger the step...a kind of logarithmic
ratio)

a=0.0001:0.0001:20;
i=(1:length(a));
x = i+i.^2/2;
x_norm = x/x(end);

now the questions:
- how do i sample then my vectors with this logarithmic sampling?
simply in this way?
K = exp(-x_norm).*sin(x_norm);

- since what i'm doing it's a kind of variable change in the integral,
do i have to take in account the Jacobian when i write
trapz(x_norm,K)

thanks.
neurino

Roger Stafford

2005-12-09, 7:27 pm

In article <1134138460.557376.16250@g43g2000cwa.googlegroups.com>,
"neurino" <lelli@physik.fu-berlin.de> wrote:

> Hallo everybody.
> I'm doomed by doubts.
>
> Let's assume i need to solve an integral of product of a vector K.
> The integrand oscillates dramatically at far right end.
> I'm *not* interested in these high frequencies, then i sample K with an
> unevenly spaced points
> (that is more far to the right, bigger the step...a kind of logarithmic
> ratio)
>
> a=0.0001:0.0001:20;
> i=(1:length(a));
> x = i+i.^2/2;
> x_norm = x/x(end);
>
> now the questions:
> - how do i sample then my vectors with this logarithmic sampling?
> simply in this way?
> K = exp(-x_norm).*sin(x_norm);
>
> - since what i'm doing it's a kind of variable change in the integral,
> do i have to take in account the Jacobian when i write
> trapz(x_norm,K)
>
> thanks.
> neurino
--------------------
Hi,

The trapz(X,Y) function takes care of non-uniformly spaced values in X.
At each interval, (x1,y1) to (x2,y2), it calculates the trapezoid area,
(y1+y2)/2*(x2-x1), which takes that into account. That is, in your
example it would give you the approximate integral of exp(-t).*sin(t) with
respect to t from t = 5e-15 to t = 1 (where t is x_norm.)

However, you worry me a little when you say, "The integrand oscillates
dramatically at far right end. I'm *not* interested in these high
frequencies, then i sample K with an unevenly spaced points ...". You
apparently have in mind going out as far as, say, t = 20 when you say
that. If it were not for the exponential decay factor, exp(-x_norm), it
would be precisely at places where the integrand "oscillates dramatically"
that you would want closer, not more widely-spaced, points to maintain
accuracy in the integration process.

(Remove "xyzzy" and ".invalid" to send me email.)
Roger Stafford

neurino

2005-12-09, 7:27 pm

However, you worry me a little when you say, "The integrand oscillates
dramatically at far right end. I'm *not* interested in these high
frequencies, then i sample K with an unevenly spaced points ...". You
apparently have in mind going out as far as, say, t = 20 when you say
that. If it were not for the exponential decay factor, exp(-x_norm),
it
would be precisely at places where the integrand "oscillates
dramatically"
that you would want closer, not more widely-spaced, points to maintain
accuracy in the integration process.
-----------------------------

Thanks for the answer, Roger.
My example above was cut just to be able to state my doubts, but the
real situation is quite different. The integrand is the product of K
and an unknown function f, that i want to retrieve, solving a linear
system.

value = trapz(a,K.*f);
min ||value - value(measured)|| in norm 2

I have some clues of its shape and regularity. It's a lognormal
distribution(PSD or SPSD). These values were taken from a model, for
testing the routines.

a = 0.0001:0.0001:20;
mu_r = [0.0690 0.3750 4.3000];
sig = [0.5641 0.6908 0.6701];

y(3,:) = N(3) .* (a(1,:) .* lognpdf(a(1,:),log(mu_r(3)),sig(3)));
y(2,:) = N(2) .* (a(1,:) .* lognpdf(a(1,:),log(mu_r(2)),sig(2)));
y(1,:) = N(1) .* (a(1,:) .* lognpdf(a(1,:),log(mu_r(1)),sig(1)));
PSD = 2 .* (sum(y));
SPSD = 4*pi .* (a.^2) .* PSD;

while the K is the sum of Bessel polynomials that converges at 2 with a
-> inf and oscillates strongly with a>3. These frequencies i guess they
are completely lost anyway in the physical process that generates my
values and, nevertheless, that minimization strongly amplify each small
perturbation in the unknown funtction.

I found out that, in order to have stability in the
minimization(inversion) above(a>4 or 5)i should have more than 100
linear systems to be solved contemporary.
And this is physically impossible for my real setup case.
Then i was approaching some new issues for me, like collocation method,
or projection method, in order to achieve the best discretization step
versus the best regularization of the linear system.
I don't still know if an uneven grid can help to stabilize my solution.
What do you think? Could i be on a false track?

neurino

Sponsored Links