Using quadgk inside a loop

How do I get quadgk to work inside a "for" loop?

I'm trying to solve an integral equation coupled with a first order ODE numerically. Therefore I need to evaluate an integral for many different steps of solving the ODE (I'm using the Euler method).

I can't seem to get my code to loop through quadgk which is the problem.

[Richard Brown]

You need to be a bit more specific about your problem - how about posting what you've tried. Then we can help.

Answer 1

[Mike Hosea]

You've got a couple of problems here that jump out at me. One is that your expression for n0 tends to lose precision due to the squaring and subtraction and cancel. Maybe rewrite it as something like n + I/2*(1-sqrt(1-4*n/I)). The second problem is that the integrand decays so rapidly that the quadrature rule does a lot of work trying to resolve it close to the left-end and then wastes a lot of energy on adding up zeros over most of the transformed interval. I suggest breaking this into a couple of pieces:

Q = quadgk(F,3.29e15,1e17,'reltol',1e-6,'abstol',1e-10) + ...
    quadgk(F,1e17,inf,'reltol',1e-6,'abstol',1e-20);

The cutoff is fairly arbitrary. It shouldn't be too large. Note here that you probably also want to crank down the absolute tolerance just a little bit on the tail because there's not much area out there to add up. -- Mike

Answer 2

[Mike Hosea]

I agree with Richard. We'll probably need to see one of your attempts, but I can take a guess. Most likely you're having trouble with the vectorization of the integrand, and it's probably erroring in mtimes something, telling you about a dimension mismatch. If that's the problem, use ARRAYFUN, e.g. if the integrand is fun(x) but fun isn't vectorized (i.e., needs to be given scalars and will not operate elementwise on an input vector), try

quadgk(@(x)arrayfun(fun,x),a,b)

This works with the new INTEGRAL function as well, and it is probably the way to go, but with INTEGRAL there's another possibility:

integral(fun,a,b,'ArrayValued',true)

Although you might not think of the integrand as returning an array, when the 'ArrayValued' option is true, the integrand is only called with scalar input.