MATLAB R2012b - Passing arrays and ints into evalin(symengine,expression)

Question

I'm attempting to generate the Laguerre polynomials, and then evaluate them elementwise over a coordinate array.

Presently my code looks something like:

[X,Y] = meshgrid(x_min:mesh_size:x_max,y_min:mesh_size:y_max);

const_p=0; 

const_l=1; %At present these two values don't really matter, any integer will do

coord_r = sqrt(X.^2 + Y.^2)

lag_input = num2str(coord_r.^2)

u_pl = evalin(symengine,['orthpoly::laguerre(',num2str(const_p),',',num2str(const_l),',',lag_input,')']);

However, that returns the following error for the last line;

Error using horzcat

Dimensions of matrices being concatenated are not consistent.

I assumed that this was because the three objects being converted to strings had different sizes, but after making them the same size the problem persists.

I'd rather avoid looping through each element if I can avoid it.

Answer 1

I would go about this slightly differently. How about the below? Note that I changed const_p and const_l from your choices because the resulting Laguerre Polynomial is spectacularly dull otherwise.

const_p = 2;
const_l = 1;

%generate the symbolic polynomial in x
lagpoly=feval(symengine,'orthpoly::laguerre',const_p,const_l,'x');

%Find the polynomical coefficients so we can evaluate using MATLAB's poly
coeff=double(feval(symengine,'coeff',lagpoly));

%generate a matrix the same size as coord_r in the original question
x=rand(512);

%Do the evaluation
u_pl=polyval(coeff,x);

Answer 2

@WalkingRandomly has the best way to do this if you need fast numeric results, which is usually the case. However, if you need exact analytical values, there is a trick that at you can use to avoid a for loop: MuPAD's map function. This is how almost all MuPAD functions must be vectorized as they're usually designed for scalar symbolic variables rather than arrays of numeric values. Here's a basic example:

const_p = 2;
const_l = 1;

mesh_size = 0.2;
x_min = 0;
x_max = 1;
y_min = 0;
y_max = 1;
[X,Y] = meshgrid(x_min:mesh_size:x_max,y_min:mesh_size:y_max);
coord_r = sqrt(X.^2 + Y.^2);

lagpoly = evalin(symengine,['map(' char(sym(coord_r)) ...
                            ',x->orthpoly::laguerre(' char(sym(const_p)) ...
                            ',' char(sym(const_l)) ',x))'])

which returns

lagpoly =

[      3,                       121/50,                   47/25,                        69/50,                   23/25,                     1/2]
[ 121/50,        76/25 - (3*2^(1/2))/5,   31/10 - (3*5^(1/2))/5, 16/5 - (3*2^(1/2)*5^(1/2))/5, 167/50 - (3*17^(1/2))/5,  88/25 - (3*26^(1/2))/5]
[  47/25,        31/10 - (3*5^(1/2))/5,   79/25 - (6*2^(1/2))/5,      163/50 - (3*13^(1/2))/5,    17/5 - (6*5^(1/2))/5, 179/50 - (3*29^(1/2))/5]
[  69/50, 16/5 - (3*2^(1/2)*5^(1/2))/5, 163/50 - (3*13^(1/2))/5,        84/25 - (9*2^(1/2))/5,                     1/2,  92/25 - (3*34^(1/2))/5]
[  23/25,      167/50 - (3*17^(1/2))/5,    17/5 - (6*5^(1/2))/5,                          1/2,  91/25 - (12*2^(1/2))/5, 191/50 - (3*41^(1/2))/5]
[    1/2,       88/25 - (3*26^(1/2))/5, 179/50 - (3*29^(1/2))/5,       92/25 - (3*34^(1/2))/5, 191/50 - (3*41^(1/2))/5,           4 - 3*2^(1/2)]

Calling double(lagpoly) will convert the result to floating point and you'll see that this is the same as the solution provided by @@WalkingRandomly (given the same inputs). Of course you could probably use the symbolic polynomial or its coefficients to find the same thing manually, though it's unfortunate that polyval isn't overloaded for class sym (there's evalp but it's also not vectorized so it would need to be used in conjunction with map).