2

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.

Community
  • 1
  • 1
Daniel Blay
  • 145
  • 7
  • Are you sure that `lag_input` isn't actually a column vector? Also, you might find `feval` more convenient. – horchler Jan 08 '14 at 04:53
  • @horchler - coord_r is a square array (presently 512x512) and lag_input has size 512x6650. I'll check out feval, cheers. – Daniel Blay Jan 08 '14 at 04:56
  • The original data size doesn't matter. You're concatenating everything horizontally in a string so the strings must all be row vectors. Check their sizes. `feval` may work if the MuPAD function has been designed for vectorized inputs - most aren't so other tricks are needed. – horchler Jan 08 '14 at 05:02
  • Each string row vector is 1x6650. I do have the suspicion that the orthpoly:laguerre function doesn't support vectorized inputs. I'm not sure what the best solution is other than to loop through the array (which seems like the least efficient method) – Daniel Blay Jan 08 '14 at 05:53

2 Answers2

2

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);
WalkingRandomly
  • 4,537
  • 5
  • 34
  • 34
  • This seems to work but it generates some weird noise (I assume from the rand), should be able to fix it. It's certainly fast! – Daniel Blay Jan 09 '14 at 01:13
  • Yes, you don't need rand at all. I just used that because I wanted a test input that was the same size as your coord_r. For your production code, remove the line `x=rand(512)` and the evaluation should look like `u_pl=polyvat(coeff,coord_r)` – WalkingRandomly Jan 09 '14 at 14:19
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).

horchler
  • 18,384
  • 4
  • 37
  • 73
  • This is a robust solution but it took me nearly an hour to run (because I'm dealing with large arrays), so I upvoted but @WalkingRandomly's solution is the ideal one for me. – Daniel Blay Jan 09 '14 at 01:13
  • @WalkingRandomly: You might also be interested in [`zip`](http://www.mathworks.com/help/symbolic/mupad_ref/zip.html) which is useful for vectorizing a function with two equal size arrays. `symobj::mapcatch` is also useful. See [here](https://github.com/horchler/SHCTools/blob/master/SHCTools/stoneholmes/private/gammaincq.m) for some examples of their usage. ;-) – horchler Jan 09 '14 at 16:49