Using quadgk inside a loop

1 次查看（过去 30 天）
David
David 2012-4-15
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.
  1 个评论
Richard Brown
Richard Brown 2012-4-16
You need to be a bit more specific about your problem - how about posting what you've tried. Then we can help.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Mike Hosea
Mike Hosea 2012-4-16
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

更多回答（1 个）

Mike Hosea
Mike Hosea 2012-4-16
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.

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by