Integral of complicated function

Hi everyone, I'm quite new to matlab, so sorry if my question may be appear silly. I want to compute numerically this integral:

where

And I want to compute this pi(a) for a = 18,19,20,...,65. Is it possible to do it? I was thinking to construct a loop and compute it for each value of a, but I have no clue on how to do it also because One of my problems is also the s variable present in the definite integral inside the definition of lambda. (gamma, c(s) are known)

Maybe it is not that difficult, but I'm totally new to matlab.

Each kind of suggestion is welcome! Thank you in advance!

[Changjie Guan]

Hi Martina, what is c(s) in lambda(s), cos(s)? To calculate this integral, I recommend you to use arrayfun().

[Martina]

no, it's a function of s.. I was lazy to plug its expression but it's a simple exponential function! And thank you for your suggestion!

Answer 1

[jgg]

There are a bunch of techniques you could use. Matlab has a built-in quadrature integration. You call it like:

 q = integral(fun,xmin,xmax)

Another way to do this is to use integration by simulation. You can look up references for the technique, but essentially you just uniformly sample R points you need to integrate. The expectation of f of these values is exactly the integral. However, this can be computationally kind of tedious.

The main issue is going to be reducing the number of integrals you want to be computing. You can break your problem apart a little bit. My workflow would be to first, just code everything up as-is to see how slow the integral is. That is, write a function for \lambda(s), c(s), and the first part of equation(1), then calculate \pi(a).

Afterwards, you can try breaking it into pieces by pre-evaluating certain integrals for values in the range you're interested in, and using trapz instead.

[jgg]

I've also noticed you can simplify lambda(s) a little bit:

 lambda(s) = c(s)*exp(s)*int_18^65(c(r)*exp(-r)dr) %not matlab code

You could write a little function

 fl = @(r)(c(r)*exp(-r));

Then, the latter part is just a number

 cons = integral(fl,18,65)

And then

lambda(s) = c(s)*exp(s)*cons;

Which will save you one integral to compute.