Vectorize MATLAB Matrices in \mathbf{v} for Spherical Heat Equation

[hunt_mat]

TL;DR

Computing a Greens function in MATLAB.

Suppose I have a vector of matrices:
<br /> \mathbf{v}=(A_{1},\cdots,A_{n})<br />

How would I vectorise this in MATLAB?

This question comes from a requirement to compute a Greens function for the spherical heat equation. I can easily compute a single function for a single position in space, but can I do it for a vector of positions, for this I think I need a vector of Greens functions for each position.

 

[RPinPA]

If the $A$'s are all the same size then v could be a 3-dimensional matrix. That would be most likely to give you a path to vectorizing. It always takes me a lot of experimenting to get things right with multi-dimensional matrices in Matlab but it's pretty sweet when it works.

Otherwise v might have to be a cell array and I'm not sure what operations might be available to do what you want to do.

Does ARRAYFUN have what you need? It seems to handle a lot of different kinds of input objects.

 

[hunt_mat]

So what I want is to compute the Greens function:
<br /> G(r,\bar{r},t)=\frac{3\bar{r}^{2}}{R^{3}}+\frac{2\bar{r}}{rR}\sum_{n=1}^{\infty}(1+\mu_{n}^{2})\sin\left(\frac{\mu_{n}\bar{r}}{R}\right)\sin\left(\frac{\mu_{n}r}{R}\right)e^{-\frac{D\mu_{n}^{2}t}{R^{2}}}<br />

From there I want to compute:

<br /> \int_{0}^{t}G(r,R,t-\tau)f(\tau)d\tau<br />

I could use a double loop no problem, but those tend to be rather slow and I wanted to vectorise it.

Mat

 

[RPinPA]

So your $A$'s are slices of $G(r,\bar r,t)$ at different $t$ values? I think this is vectorizable with judicious use of MESHGRID to generate all the combinations of the input arguments, SUM over particular dimensions and CONVN to generate the convolution at the end.

As I said, I always have to do a lot of experimentation to get multidimensional calculations like this right. For instance what order to I want the inputs to MESHGRID? But I'm pretty sure those three functions have what you need to do this.

I will in fact do a little experimenting later for more specific guidance.

 

[hunt_mat]

No, the A's are slices of G(r,\bar{r},t) at different r's. As I said, I can do it using a double loop, but it's slow like that, I wanted speed as I'll be using it lots and lots of times.

Mat

 

[RPinPA]

OK, I did a little experimenting as promised.

I started with this: [r, rbar, t, n] = ndgrid(r values, rbar values, t values, n values). You may want to use a different order of arguments.

You can then evaluate the expression under your summation sign using element-wise multiplication, i.e. B = sin(...r...) .* sin(...rbar...).*exp(...t...), etc.

Use SUM(B, 4) to sum over the index n, which is the 4th index for my example.

To finish evaluating the expression for G, you can access subsets of r and rbar for one value of n, i.e. r(:, :, :, 1) and rbar(:, :, :, 1).

That gives the entire 3-D array $G(r, \bar r, t)$ with $r$ being the first index, $\bar r$ being the second and $t$ being the third.

The convolution can then be done with CONVN, but I think you may need to build an array representing $f(\tau)$ which is the same shape as G. I'm a little confused from a quick reading on what exactly CONVN is doing when one argument is 3-dimensional and another is 1-dimensional.

But the bottom line is I'm 99% sure it can be vectorized in this way.

 

[hunt_mat]

Okay, cheers. I'll play with that and see if it works.

 

[RPinPA]

I'll send you the transcript of my experiments if you're interested.

 

[hunt_mat]

Please. I've been ill with a very bad cold and haven't been able think straight but I'm better (not fully) to begin thinking again.

 

[hunt_mat]

I can't do this via vectorisation can I? It has to be for loops.

 

[hunt_mat]

This is the function I have:

[CODE lang="matlab" title="Greens function"]function G = Greens_2(D,mu_n,r,t,r_bar,n)
%This is the Green's function for spherical diffusion equation
%Solution can be found in A.D. Polyanin. The Handbook of Linear Partial Differential Equations
%for Engineers and Scientists, Chapman & Hall, CRC 17
n=length(mu_n);
N=1:n;
G_n=zeros(length(t),length(r},length(r_bar),n);
R=r_bar(end);
N=length(mu_n); %Number of terms used in Green's function series
N_r=length(r_bar);
for l=1:n
for k=1:N_r
for j=1:N_r
G_n(:,j,k,l)=(2*r_bar(k)/r(j)*R)*(1+(mu_n(l))^(-2))*sin(mu_n(l)*r_bar(k)/R)*sin(mu_n(l)*r(j)/R)*exp(-D*mu_{n}(l}^2*t/R^2);
end
end
end

G=sum(G_n,4)+3*r_bar.^2/R^3;[/CODE]