求问matlab下用循环建立符号数组的办法

Question

RT，现在希望实现这样一个功能：需要创建一个较大的四维数组（50*50*50*50），每个元素都是类似sin((i+j+k+l)*pi*x)这种形式，其中x为符号，请问该用哪个指令？多谢。
顺带一提，输入完之后还要对这个数组做定积分，int命令会不会太慢？如果是的话，还请大大们提个建议的算法。

Answer 1

就这个问题自身而言，并不算困难，例如，用下面几行代码就可以：

N=5;
[i,j,k,l]=ndgrid(1:N,1:N,1:N,1:N);
syms x
y=sin((i+j+k+l)*pi*x);
I=int(y,0,0.5);

注意，这里我只是示范解决问题的基本思路，所以，把问题规模大幅度减小为N=5。

 

使用MATLAB的profiler可以分析执行这些代码所需的时间，以下是取N=5和10所对应的结果：

 

我另外尝试了N=20的情况，在我的机器上算了1个小时还没出来结果。具体情况和硬件配置、操作系统、MATLAB版本等因素都有关系，所以，这些只能作为一个参照。尽管如此，我们也可以想象，如果取N=500，即使现在最豪华配置的PC，计算时间只怕也要以天为单位（而且，很难保证是否会出现内存溢出等情况）。

 

通过上面的分析，基本上可以说，这种做法不可行。有没有可能通过对程序进行优化来解决呢？上面的代码使用了向量化的运算，照理说，已经是MATLAB最高效率的代码写法了。如果改写成显式循环的方法，代码如下：

N=5;
[i,j,k,l]=ndgrid(1:N,1:N,1:N,1:N);
syms x
y=sym(zeros(N,N,N,N));
n=i+j+k+l;
I=zeros(N,N,N,N);
for ii=1:N^4
    y(ii)=sin(n(ii)*pi*x);
    I(ii)=int(y(ii),0,0.5);
end

这段代码对于N=5和N=10的运行结果如下：

 

可见，使用循环的代码比起向量化代码，所需时间还要长很多。

 

那么，这个问题究竟还有没有更有效率的解法呢？

 

我认为，如果每个元素都是类似sin((i+j+k+l)*pi*x)这种形式，比较合理的做法应该是，把系数(i+j+k+l)看成一个常数，记作n，先对sin(n*pi*x)求定积分得到关于n的一般表达式，然后再把具体的系数代入，代码大致如下：

N=50;
syms x n
f=sin(n*pi*x);
I=int(f,0,0.5);
[i,j,k,l]=ndgrid(1:N,1:N,1:N,1:N);
I=inline(I);
I=I(i+j+k+l);

注意，最后两行先转换为inline函数再代入，而不是直接subs，二者效率差别十分巨大，可自行比较。

 

最后看一下结果。这里即使取N=50（如果N更大，所需内存较多），所需时间也很短，而前面的算法即使对于N=20已经基本上无法计算了：

上述代码在2007b上调试，在其他版本可能会有不同，特此提示。如果有疑问可继续追问。

追问

多谢多谢。这个问题其实我已经用另一个方法解决了。
虽然笨了点，但很有效。
syms x;
B=zeros(n,n,n,n);
A=sym(B);下面的循环打不下了，坑爹，但反正就是这样解决了问题。
另外，50*50*50*50的数组用matlab计算不太现实，关键是符号变量太大，把内存耗尽了，我问了问老师，他说取到20也能满足要求精度。
话说,QUADV指令效率高么？

追答

50*50*50*50的数组确实比较大，如果是double类型，需要50M的内存，但sym所需的内存要多的多，而且计算过程中要启动符号运算引擎（Maple或MuPad），相关的数据在运算引擎内也有相应的内存需求（很可能不少于MATLAB自身），只怕现在大多数的PC机还很难承受。

就你的问题自身而言，因为各元素存在通式，所以可以先用符号运算把积分的函数表达式求出来，然后使用inline函数直接计算，这样速度快而且精度有保证。可能的条件下尽量用这种方法。

quadv是基于纯数值计算的，效率应该还可以，要看具体的问题和做法。不过，可以肯定，效率肯定远远不如上面说的先符号后inline的方法，但比起直接用int计算效率会高。