Matlab - questions about numerical integrals

Author

questions about numerical integrals

Vista

2007-06-25, 10:15 pm

Hi all,

Recently I was puzzled a lot by the numerical integrals.

1. When doing the inverse Laplace Transform via the Bromwich integral,

integrate( f(s)*exp(-s*x), s from c-i*inf to c+i*inf),

theoretically, c doesn't matter as long as c is on the right side of all the
singularities.

But numerically, c affects the numerical results and it varies a lot.

And since I don't have a closed-form result to compare with the numerical
results, I just couldn't tell if the numerical results are correct or
not. Even worse, I couldn't tell which "c" gives me the correct value, as
the results differ a lot depending on different values of "c". If I plot
the value of numerical integral vs. the value of "c", there are spikes,
i.e. for a certain "c", the numerical integral gives huge spike/peak with no
reasons.

Does anybody know how to find an "optimal" "c"?

I've attached my function at the end of this post, please take a look! Thank
you very much for your help!

2. When computing the integral from -inf to +inf, of course I use large
values to truncate the integral, for example, I set L=100000 and did
integrate( f(c+i*u)*exp(-(c+i*u)*x), u from -L to L)
In
quadl(@(u) f(c+1i*u)*exp(-(c+1i*u)*x)*1i, -L, L)

But such truncation does not stabalize for large L's. Sometimes large L's
even give values that are way off, than small L's. And large L's tend to
give explosive numerical integral results, even for f(.) decays to zero very
rapidly.

3. When doing the transformation of integral interval from [-inf, +inf]
to [-1, 1] using y=tanh(x), for most values of "c", the transformed
finite interval numerical integral even gives worse result than the
nontransformed infinite interval numerical integral(the original one),
except a small number of values of "c". Of course for this particular set
of parameter combination, I can decide the optimal value of "c", which
gives me smaller error for the transformed integral than the
non-transformed integral. But if I change the parameters, everything
changed and there seems to be no "automatic" way to handle this problem.

4. Here is a result from Matlab:

>Warning: Maximum function count exceeded; singularity likely.
> In quadl at 98

What was the problem?

Is there a way I can find out where is the singularity point? Because it
is very complicated, so when I did visual inspection at theory level, I
didn't find out any singularity, also due to its complex-valued natures.

5. I desperately need reference on oscillatory integrals and infinite
interval
integrals. The pain of dealing with numerical integrations is that you don't
have closed form or theoretically correct result to compare to. Either you
trust the numerical results, or you find them fishy but you don't know
what's wrong...

Thanks a lot for your help!

p.s. For this piece of code, change-of-variable of y=tanh(x) seems works,
but I don't know how to stress test its correctness for all parameters
combinations.

Yet, still, how to choose the "optimal" "c"?

--------------------------------
t=6;
y0=1;
theta=1;
kappa=1;
lambda=2;
eta=[2; 4; 7];
rho=[1/3; 1/3; 1/3];
L = 1000000;
c= 0.1;
MyIntResult=quadl(@(x) real(myIntegrand(x, t, theta, y0, kappa, lambda, eta,
rho, c, s)), -L, L)/2/pi;

--------------------------------

function res=myIntegrand(xs, t, theta, y0, kappa, lamda, eta, rho, C, s)

%rho and eta are vectors of same length.

for i=1:length(xs);
x=xs(i);
u=-(C+1i*x);
res(i)=exp(-u*theta*t-u*(y0-theta)/kappa*(1-exp(-kappa*t))-
rho'*(lamda*u*t./(u+kappa*eta)) ...
+rho'*(lamda./(u+kappa*eta).*log(1+u./eta/kappa*(1-exp(-kappa*t)))));
res(i)=res(i)*exp(s*u)/(-u);
end;

Michael Jørgensen

2007-06-26, 8:13 am

"Vista" <abc@gmai.com> wrote in message
news:f5pveq$1dm$1@news.Stanford.EDU...
> Hi all,
>
> Recently I was puzzled a lot by the numerical integrals.
>
> 1. When doing the inverse Laplace Transform via the Bromwich integral,
>
> integrate( f(s)*exp(-s*x), s from c-i*inf to c+i*inf),
>
> theoretically, c doesn't matter as long as c is on the right side of all
> the singularities.
>
> But numerically, c affects the numerical results and it varies a lot.

1) Well, if your numerical integration gives different results depending on
the value of c, then your numerical integration is incorrect OR c is not on
the right of all the singularities.

2) Integrating over a very large interval, like +- 1e6, it is very likely
the numerical integration will miss out any "interesting details" on a much
smaller scale. I would instead calculate the integral of, say, [-10,10].
Then I would calculate [10,20] and see if that value is much smaller than
the main part. If not, I'll keep calculating [10 k, 10 (k+1)], for integer
k, and when k increases, the integral over this interval should converge
rapidly towards zero.

3) The substitution y = tanh(x) is a good idea to get rid of the infinite
interval. However, there may still be singularities at the end points y = +-
1. So you should find out, if there is a way to "help" the numerical
integration around these singularities. Or at least to increase the accuracy
of the numerical integration.

-Michael.

toni.lassila@gmail.com

2007-06-26, 8:13 am

On Jun 26, 6:00 am, "Vista" <a...@gmai.com> wrote:
> Hi all,
>
> Recently I was puzzled a lot by the numerical integrals.
>
> 1. When doing the inverse Laplace Transform via the Bromwich integral,
>
> integrate( f(s)*exp(-s*x), s from c-i*inf to c+i*inf),
>
> theoretically, c doesn't matter as long as c is on the right side of all the
> singularities.
>
> But numerically, c affects the numerical results and it varies a lot.

Unfortunately, the numerical computation of the inverse Laplace
transform is well known to be an ill-posed problem so that any simple
numerical method you try will probably run into some problems.

Rather than trying to hack together your own method, why not study
some of the many papers written on the topic? You can start from:

http://www.google.com/search?hl=en&...sform
%22

widmar

2007-06-29, 10:12 pm

Vista wrote:
>
> 5. I desperately need reference on oscillatory integrals and infinite interval
> integrals. The pain of dealing with numerical integrations is that you don't
> have closed form or theoretically correct result to compare to. Either you
> trust the numerical results, or you find them fishy but you don't know
> what's wrong...

quadpack (netlib) has a sizable collection of quadrature codes, so you
should start looking there. To get you started,

- qagi.f (for infinite intervals)
- qag.f (for oscillating integrands)

Regarding pain, it could be largely due to your using a blackbox sw -
it's not too hard to find or construct a case with known analytical
properties which can then be used to "stress" test a given solver.
e.g. see the last post in a "Is there a good way to compute indefinite
integral..." thread.

---
sdx - modeling, simulation.
http://www.sdynamix.com

Sponsored Links