Thread Subject:

quad error

Hello,

I have a simple function to numerically integrate and it depends on
the scalar 't':

%%%%%%%%%%%%%%%%%
function omegafun = omegafun(mM1, t)
 omegafun = sqrt(det(eye(size(mM1, 1)) -2*t.*mM1))
end
%%%%%%%%%%%%%%%%%

Here mM1 is a parameter matrix. I am able to compute this function and
plot it, but when I try quad,

quad(@(t) omegafun(mM1, t), 0, 10);

I get the error:

??? Error using ==> times
Matrix dimensions must agree.

Error in ==> omegafun at 2
 omegafun = sqrt(det(eye(size(mM1, 1)) -2*t.*mM1))

Error in ==> @(t)omegafun(t,mM1)


Error in ==> quad at 77
y = f(x, varargin{:});


Something obvious that I am missing?

"tzygmund" <tzygmund@googlemail.com> wrote in message
news:2d350472-226f-4122-8c33-4ed7b7ee5ddc@t19g2000prd.googlegroups.com...
>
> Hello,
>
> I have a simple function to numerically integrate and it depends on
> the scalar 't':
>
> %%%%%%%%%%%%%%%%%
> function omegafun = omegafun(mM1, t)
> omegafun = sqrt(det(eye(size(mM1, 1)) -2*t.*mM1))
> end
> %%%%%%%%%%%%%%%%%
>
> Here mM1 is a parameter matrix. I am able to compute this function and
> plot it, but when I try quad,
>
> quad(@(t) omegafun(mM1, t), 0, 10);

http://www.mathworks.com/help/techdoc/ref/quad.html

"The function y = fun(x) should accept a vector argument x and return a
vector result y, the integrand evaluated at each element of x."

From your description it sounds like your function cannot handle a vector
input t. Modify it so that it can by doing something like looping over the
elements of t. Don't forget to preallocate omegafun to be the same size and
shape as t.

omegafun = zeros(size(t));

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

tzygmund <tzygmund@googlemail.com> wrote in message <2d350472-226f-4122-8c33-4ed7b7ee5ddc@t19g2000prd.googlegroups.com>...
> ........
> I have a simple function to numerically integrate and it depends on
> the scalar 't':
>
> %%%%%%%%%%%%%%%%%
> function omegafun = omegafun(mM1, t)
> omegafun = sqrt(det(eye(size(mM1, 1)) -2*t.*mM1))
> end
> %%%%%%%%%%%%%%%%%
> .......
- - - - - - - - - -
  Your integrand is the square root of a polynomial in t that can be expressed in terms of the eigenvalues of nM1. You can use this fact to allow t to be accepted as a vector.

[D,T] = ndgrid(eig(nM1),t);
omega = sqrt(prod(1-2*T.*D,1));

The above will accept the argument t as a vector. (Write your function so it accepts the eigenvalues of nM1 as input rather than nM1 itself to keep from repeatedly computing the same eigenvalues.)

  (It is to be hoped that the range of t will keep the above product always non-negative to avoid imaginary values for the square root. I am not sure 'quad' would handle them properly.)

Roger Stafford

Thanks Roger and Steven for excellent suggestions.


On Apr 19, 4:18 pm, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> tzygmund <tzygm...@googlemail.com> wrote in message <2d350472-226f-4122-8c33-4ed7b7ee5...@t19g2000prd.googlegroups.com>...
> > ........
> > I have a simple function to numerically integrate and it depends on
> > the scalar 't':
>
> > %%%%%%%%%%%%%%%%%
> > function omegafun = omegafun(mM1, t)
> >  omegafun =  sqrt(det(eye(size(mM1, 1)) -2*t.*mM1))
> > end
> > %%%%%%%%%%%%%%%%%
> > .......
>
> - - - - - - - - - -
>   Your integrand is the square root of a polynomial in t that can be expressed in terms of the eigenvalues of nM1.  You can use this fact to allow t to be accepted as a vector.
>
> [D,T] = ndgrid(eig(nM1),t);
> omega = sqrt(prod(1-2*T.*D,1));
>
> The above will accept the argument t as a vector.  (Write your function so it accepts the eigenvalues of nM1 as input rather than nM1 itself to keep from repeatedly computing the same eigenvalues.)
>
>   (It is to be hoped that the range of t will keep the above product always non-negative to avoid imaginary values for the square root.  I am not sure 'quad' would handle them properly.)
>
> Roger Stafford