Integrals in Freemat

[george]

Hello,

How do I take the integral of a function in Freemat? I know Matlab
has
the int() function.

Thank you.

[Eugene]

FreeMat does not support symbolic operations (yet). You need symbolic
algebra software, e.g. Maxima (http://maxima.sourceforge.net/).

[Timothy Cyders]

If a numerical solution will suit you (indefinite integrals need not apply), you might look at the cumtrapz() function. Example: if you want to take the integral of y = x^2+4*x from 0 to 10,

x = 0:0.01:10; % better resolution will give you more precise results
y = x.^2+4*x;
integral = cumtrapz(x,y) % gives a vector of values for the definite integral from x(1) to x(n). The last value integral(length(integral)) is the overall definite integral.

TJ

On Tue, Nov 16, 2010 at 2:38 AM, Eugene <gen...@gmail.com> wrote:

FreeMat does not support symbolic operations (yet). You need symbolic
algebra software, e.g. Maxima (http://maxima.sourceforge.net/).

--
You received this message because you are subscribed to the Google Groups "freemat" group.
To post to this group, send email to fre...@googlegroups.com.
To unsubscribe from this group, send email to freemat+u...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/freemat?hl=en.

[jonw0224]

I don't have the cumtrapz function with FreeMat. I looked at Octave
and it's command isn't the same as MATLAB's. I wrote my own that
seems to work up to three dimensions anyway. Here it is inline:

% CUMTRAPZ - Cumulative trapezoidal numerical integration
%
% Usage
%
% Z = cumtrapz(Y)
% Z = cumtrapz(X,Y)
% Z = cumtrapz(Y,dim)
% Z = cumtrapz(X,Y,dim)
%
% Description
%
% Z = cumtrapz(Y) computes an approximation of the cumulative integral
of Y via
% the trapezoidal method with unit spacing. To compute the integral
with other
% than unit spacing, multiply Z by the spacing increment. Input Y can
be complex.
%
% For vectors, cumtrapz(Y) is a vector containing the cumulative
integral of Y.
%
% For matrices, cumtrapz(Y) is a matrix the same size as Y with the
cumulative
% integral over each column.
%
% For multidimensional arrays, cumtrapz(Y) works across the first
% nonsingleton dimension.
%
% Z = cumtrapz(X,Y) computes the cumulative integral of Y with respect
to X
% using trapezoidal integration. X and Y must be vectors of the same
length,
% or X must be a column vector and Y an array whose first nonsingleton
dimension
% is length(X). cumtrapz operates across this dimension. Inputs X and
Y can be
% complex.
%
% If X is a column vector and Y an array whose first nonsingleton
dimension
% is length(X), cumtrapz(X,Y) operates across this dimension.
%
% Z = cumtrapz(X,Y,dim) or cumtrapz(Y,dim) integrates across the
dimension of
% Y specified by scalar dim. The length of X must be the same as
size(Y,dim).
%
% Licensed under the GPL
% Copyright (C) 2010 Jonathan Weaver <jonw...@netscape.net>

function [A] = cumtrapz(x,y,dim)

% This program is free software; you can redistribute it and/or modify
it under
% the terms of the GNU General Public License as published by the Free
Software
% Foundation; either version 2 of the License, or (at your option) any
later
% version.
%
% This program is distributed in the hope that it will be useful, but
WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for
more
% details.
%
% You should have received a copy of the GNU General Public License
along with
% this program; if not, write to the Free Software Foundation, Inc.,
% 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

if (nargin == 1) %cumtrapz(y)
[y, dim] = shiftdim(x);
x = 1:size(y,1);
elseif(nargin == 2) %cumtrapz(x,y) or cumtrapz(y,dim)
if(isscalar(y)) %cumtrapz(y,dim)
dim = y - 1;
y = shiftdim(x, dim);
x = 1:size(y,1);
else %cumtrapz(x,y);
[y, dim] = shiftdim(y);
if(length(x) ~= size(y,1))
error('The length of x must be equal to the length of
the first non-singleton dimension of y');
end
end
elseif(nargin == 3) %cumtrapz(x,y,dim)
dim = dim - 1;
y = shiftdim(y,dim);
if(length(x) ~= size(y,1))
error('The length of x must be equal to the size(Y,dim)');
end
else
error('Invalid number of arguments');
end
n = size(y);
m = n(2:end);
a = [zeros([1 m]); (y(1:(n(1)-1),:,:,:,:,:)+y(2:n(1),:,:,:,:,:))/
2.*repmat(diff(x)', [1 m])];
A = shiftdim(cumsum(a),length(n) - dim);
end

On Nov 16, 6:02 am, Timothy Cyders <t.cyd...@gmail.com> wrote:
> If a numerical solution will suit you (indefinite integrals need not apply),
> you might look at the cumtrapz() function. Example: if you want to take the
> integral of y = x^2+4*x from 0 to 10,
>
> x = 0:0.01:10; % better resolution will give you more precise results
> y = x.^2+4*x;
> integral = cumtrapz(x,y) % gives a vector of values for the definite
> integral from x(1) to x(n). The last value integral(length(integral)) is the
> overall definite integral.
>
> TJ
>

> On Tue, Nov 16, 2010 at 2:38 AM, Eugene <gene...@gmail.com> wrote:
> > FreeMat does not support symbolic operations (yet). You need symbolic
> > algebra software, e.g. Maxima (http://maxima.sourceforge.net/).
>
> > --
> > You received this message because you are subscribed to the Google Groups
> > "freemat" group.
> > To post to this group, send email to fre...@googlegroups.com.
> > To unsubscribe from this group, send email to

> > freemat+u...@googlegroups.com<freemat%2Bunsu...@googlegroups.com>
> > .

[Timothy Cyders]

Cumtrapz is included in the latest version of FreeMat - I wrote it into the code myself. What version are you using?

TJ

To unsubscribe from this group, send email to freemat+u...@googlegroups.com.

[jonw0224]

TJ,

It looks like I hadn't downloaded the latest version. I'm now using
FreeMat v4.0 for windows. I see the cumtrapz.m file in the matrix
directory.

First of all, thank you for working on such a great piece of
software. However, cumtrapz that's included with v4.0 isn't
functionally the same as MATLAB's description at
http://www.mathworks.com/help/techdoc/ref/cumtrapz.html, which I was
attempting to duplicate. Also, performance is lacking due to the
interpreted loops. I'd encourage you to try my implementation and I'm
more than happy to contribute it to the project.

For performance testing, I renamed my implementation as cumtrapz2.m
and typed the following in an m file:

x = 0:0.01:100;
y = sin(x);
tic;
z = cumtrapz(x,y);
a = toc
tic;i
z2 = cumtrapz2(x,y);
b = toc
issame(z,z2) % I'm having trouble with isequal in the new version,
worked in the other version???

Executing it, I get

a =
66.2980
b =
0.1250
ans =
1

I get much worse results when I change the first line to x =
0:0.01:1000 ...

-Jonathan

[Timothy Cyders]

Jonathan,

 

Your function looks good. If I recall correctly, mine calls a loop in such a way that it cannot be compiled by JIT. I actually wrote that one mostly on paper in Africa. Yours looks much better without all the loops and with matrix operations instead. I'll do some testing on it and commit it to the code as cumtrapz() when I get the chance. Thanks for the contribution!

 

TJ

--

[jonw0224]

TJ,

I just can't leave code alone. I thought that I'd try to write the
trapz function since it is similar to cumtrapz. While writing the
trapz function, I kept running into problems trying to visualize 3
dimensional matrices which led to many bugs. It hit me last night as
I was going to sleep to use reshape to flatten the matrix to 2
dimensions and then I applied that idea to the implementation
cumtrapz. It actually improves the code readability a good bit and
performance slightly. On a related note, is this group the best way
to share/contribure m-files written for FreeMat? For now, both
functions are inline below.

-Jonathan

[y, dim] = shiftdim(x); %rotate matrix

x = 1:size(y,1);
elseif(nargin == 2) %cumtrapz(x,y) or cumtrapz(y,dim)
if(isscalar(y)) %cumtrapz(y,dim)
dim = y - 1;

y = shiftdim(x, dim); %rotate matrix

x = 1:size(y,1);
else %cumtrapz(x,y);

[y, dim] = shiftdim(y); %rotate matrix

if(length(x) ~= size(y,1))
error('The length of x must be equal to the length of
the first non-singleton dimension of y');
end
end
elseif(nargin == 3) %cumtrapz(x,y,dim)
dim = dim - 1;

y = shiftdim(y,dim); %rotate matrix

if(length(x) ~= size(y,1))
error('The length of x must be equal to the size(Y,dim)');
end
else
error('Invalid number of arguments');
end
n = size(y);

m = n(1);
p = prod(n(2:end));
y = reshape(y, [m p]); % flatten matrix
b = (y(1:(m-1),:) + y(2:m,:)) / 2;
c = repmat(diff(x).', [1 p]);
a = [zeros([1 p]); cumsum(b.*c)];
a = reshape(a, n); % unflatten matrix
A = shiftdim(a,length(n) - dim); % unrotate matrix
end


% TRAPZ - Trapezoidal numerical integration
%
% Usage
%
% Z = trapz(Y)
% Z = trapz(X,Y)
% Z = trapz(Y,dim)
% Z = trapz(X,Y,dim)
%
% Description
%
% Z = trapz(Y) computes an approximation of the cumulative integral of

Y via
% the trapezoidal method with unit spacing. To compute the integral
with other
% than unit spacing, multiply Z by the spacing increment. Input Y can
be complex.
%

% For vectors, trapz(Y) is the integral of Y.
%
% For matrices, cumtrapz(Y) is a row vector with the integral over
each column.
%
% For multidimensional arrays, trapz(Y) works across the first
% nonsingleton dimension.
%
% Z = trapz(X,Y) computes the cumulative integral of Y with respect to

X
% using trapezoidal integration. X and Y must be vectors of the same
length,
% or X must be a column vector and Y an array whose first nonsingleton
dimension

% is length(X). trapz operates across this dimension. Inputs X and Y

can be
% complex.
%
% If X is a column vector and Y an array whose first nonsingleton
dimension
% is length(X), cumtrapz(X,Y) operates across this dimension.
%

% Z = trapz(X,Y,dim) or trapz(Y,dim) integrates across the dimension

of
% Y specified by scalar dim. The length of X must be the same as
size(Y,dim).
%
% Licensed under the GPL
% Copyright (C) 2010 Jonathan Weaver <jonw...@netscape.net>

function [A] = trapz(x,y,dim)

% This program is free software; you can redistribute it and/or modify
it under
% the terms of the GNU General Public License as published by the Free
Software
% Foundation; either version 2 of the License, or (at your option) any
later
% version.
%
% This program is distributed in the hope that it will be useful, but
WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for
more
% details.
%
% You should have received a copy of the GNU General Public License
along with
% this program; if not, write to the Free Software Foundation, Inc.,
% 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

if (nargin == 1) %trapz(y)
[y, dim] = shiftdim(x); %rotate matrix
x = 1:size(y,1);
elseif(nargin == 2) %trapz(x,y) or trapz(y,dim)
if(isscalar(y)) %trapz(y,dim)

dim = y - 1;

y = shiftdim(x, dim); % rotate matrix
x = 1:size(y,1);
else %trapz(x,y);
[y, dim] = shiftdim(y); % rotate matrix

if(length(x) ~= size(y,1))
error('The length of x must be equal to the length of
the first non-singleton dimension of y');
end
end

elseif(nargin == 3) %trapz(x,y,dim)

dim = dim - 1;

y = shiftdim(y,dim); % rotate matrix

if(length(x) ~= size(y,1))
error('The length of x must be equal to the size(Y,dim)');
end
else
error('Invalid number of arguments');
end
n = size(y);

m = n(1);
m2 = n(2:end);
y = reshape(y, [m prod(m2)]); % flatten matrix
b = (y(1:(m-1),:) + y(2:m,:)) / 2;
c = diff(x);
a = c * b;
a = reshape(a, [1 m2]); % unflatten matrix
A = shiftdim(a,length(n) - dim); % rotate matrix
end

On Nov 18, 12:34 pm, Timothy Cyders <t.cyd...@gmail.com> wrote:
> Jonathan,
>
> Your function looks good. If I recall correctly, mine calls a loop in such a
> way that it cannot be compiled by JIT. I actually wrote that one mostly on
> paper in Africa. Yours looks much better without all the loops and with
> matrix operations instead. I'll do some testing on it and commit it to the
> code as cumtrapz() when I get the chance. Thanks for the contribution!
>
> TJ
>
> On Thu, Nov 18, 2010 at 11:36 AM, jonw0224

> <jonathan.wea...@windstream.net>wrote:

> > freemat+u...@googlegroups.com<freemat%2Bunsu...@googlegroups.com>
> > .

[Paul Picot]

On Fri, Nov 19, 2010 at 09:23, jonw0224 <jonatha...@windstream.net> wrote:
>  On a related note, is this group the best way
> to share/contribure m-files written for FreeMat?

I know it's not the best way to contribute code in a collaborative
project but I, for one, appreciate seeing discussion and sharing of
code snippets in this forum. It's a great way to share best-practices
and learn.

[Timothy Cyders]

I would tend to agree, what do you think Eugene?

TJ

--
You received this message because you are subscribed to the Google Groups "freemat" group.
To post to this group, send email to fre...@googlegroups.com.

To unsubscribe from this group, send email to freemat+u...@googlegroups.com.

[Eugene]

I like list discussions very much. However, the best way to make sure
that the patch makes it into the next release is to submit it as a bug
or a patch on the sourceforge site:
http://sourceforge.net/projects/freemat/support

Before the release Samit or I go through that list of bugs and
patches. I sometimes miss some of the discussions on the list, and
would most likely forget about a patch a few weeks later.

Thank you for contributing the code.

Eugene

On Nov 19, 6:56 am, Timothy Cyders <t.cyd...@gmail.com> wrote:
> I would tend to agree, what do you think Eugene?
>
> TJ
>
> On Fri, Nov 19, 2010 at 9:55 AM, Paul Picot <pic...@gmail.com> wrote:

> > On Fri, Nov 19, 2010 at 09:23, jonw0224 <jonathan.wea...@windstream.net>

> > wrote:
> > >  On a related note, is this group the best way
> > > to share/contribure m-files written for FreeMat?
>
> > I know it's not the best way to contribute code in a collaborative
> > project but I, for one, appreciate seeing discussion and sharing of
> > code snippets in this forum.  It's a great way to share best-practices
> > and learn.
>
> > --
> > You received this message because you are subscribed to the Google Groups
> > "freemat" group.
> > To post to this group, send email to fre...@googlegroups.com.
> > To unsubscribe from this group, send email to

> > freemat+u...@googlegroups.com<freemat%2Bunsu...@googlegroups.com>
> > .

[jonw0224]

Ok, I'll add a patch and make a posting...thanks!