一个递归或迭代求解问题，超大运算量，如何实现？

wanglu

原帖位置：http://www.ilovematlab.cn/thread-60022-1-1.html
问题描述：


我用Forcal给出了递归解法，但当m＝4时，运算时间就很长了，这倒不是Forcal速度太慢的缘故，此例若转换成C/C++代码，Forcal的速度估计在C/C++速度的1/10左右，甚至更快一些。速度慢的原因主要在于存在大量的函数调用。

如果能保存V(t,m)的计算结果，以后不再计算，估计能大幅提高速度。但难点在于，此问题中存在一个积分（我用的积分方法是徐士良算法中的连分式积分法），不知道积分运算时V(t,m)中的t为多少，故难以保存中间结果。

Matlab、mmtc、Maple等高手们看一下有什么好方法，最好给出代码，实际运行一下。

据原贴主xl0418讲，有人算过m＝50，以运行10分钟为限，大家看能算到多少？

下面是Forcal代码及结果：

1. 
   
2. SetTalk[false];
   
3. !SetRealStackMax(1000);
   
4. Rm(t,m:j,s,a,b,h:hVm,iRm)=
   
5. {
   
6.   iRm++, h=0.00001,
   
7.   sys::call(2,hVm :t+h,m : &a),  //计算s=Vm(t+h,m)
   
8.   sys::call(2,hVm :t,m : &b),    //计算s=Vm(t,m)
   
9.   s=[a-b]/h,
   
10.   j=0,(j<=m).while{
    
11.     sys::call(2,hVm :t,j : &a), sys::call(2,hVm :t,m-j : &b),  //计算a=Vm(t,j),   b=Vm(t,m-j)
    
12.     s=s+a*b,
    
13.     j++
    
14.   },
    
15.   s-[1-which(m<1,0,1)] //最后返回这个值。看这句有没有问题：取m<1，而不是m<=1
    
16. };
    
17. //Fkappa将返回s的值。mm是模块变量，因调用Rm(t,mm-1)时mm将改变，故用a预先保存，函数返回前恢复mm的值
    
18. Fkappa(t:a,s,kappa:mm,iFkappa)= iFkappa++, kappa=0, a=mm, s=t^(2*kappa)*Rm(t,mm-1), mm=a, s;
    
19. Vm(t,m:h,kappa:mm,hFkappa,iVm)=
    
20. {
    
21.   iVm++, h=-1, kappa=0,
    
22.   if[m<=0,return(t)], if[m==1,return(1/3*h*t^3)],
    
23.   which(m<=1,0,1)*Vm(t,m-1)+h*{mm=m, XSLSF::fpqg[hFkappa,0,t,1e-6]}  //XSLSF::fpqg即徐士良算法中的连分式积分法
    
24. };
    
25. pert(t)=t-1/3*t^3+2/15*t^5-17/315*t^7;
    
26. main(t,m:k,s:hFkappa,hVm,iRm,iFkappa,iVm)=
    
27. {
    
28.   iRm=0, iFkappa=0, iVm=0,
    
29.   hFkappa=HFor("Fkappa"), hVm=HFor("Vm"),  //预先获得函数Fkappa和Vm的句柄，函数调用时将用到
    
30.   s=0,k=0,(k<=m).while{
    
31.     s=s+Vm(t,k),
    
32.     k++
    
33.   },
    
34.   printff{"\r\n结果＝{1,r,-25.16}, 验证pert(t)＝{2,r,-25.16}, 函数Rm调用次数{3,i,-8}, 函数Fkappa调用次数{4,i,-8}, 函数Vm调用次数{5,i,-8}",s,pert(t),iRm,iFkappa,iVm},
    
35.   s  //返回s
    
36. };
    
37. main[1,1];
    
38. main[1,2];
    
39. main[1,3];
    
40. main[3,1];
    
41. main[3,2];
    
42. main[3,3];
    

复制代码

结果：

结果＝0.6666666666666667      , 验证pert(t)＝0.7460317460317461      , 函数Rm调用次数0      , 函数Fkappa调用次数0      , 函数Vm调用次数2      
结果＝0.8000050003347501      , 验证pert(t)＝0.7460317460317461      , 函数Rm调用次数65      , 函数Fkappa调用次数65      , 函数Vm调用次数394   
结果＝0.7460312361721437      , 验证pert(t)＝0.7460317460317461      , 函数Rm调用次数33702  , 函数Fkappa调用次数33702  , 函数Vm调用次数203761  
结果＝-6.                      , 验证pert(t)＝-91.62857142857143      , 函数Rm调用次数0      , 函数Fkappa调用次数0      , 函数Vm调用次数2      
结果＝26.40004500009          , 验证pert(t)＝-91.62857142857143      , 函数Rm调用次数129    , 函数Fkappa调用次数129    , 函数Vm调用次数778   
结果＝-91.62900091017789      , 验证pert(t)＝-91.62857142857143      , 函数Rm调用次数163158  , 函数Fkappa调用次数163158  , 函数Vm调用次数982033

我还有一个问题，放在maple板块：一个含积分的等式优化题
位置：http://forum.simwe.com/thread-910945-1-1.html
大家帮忙解决一下？

messenger

看了一下，问题描述的不是很清楚，当k=0是时，v(t,k)为v(t,0)，也就是要求第一、二两个式子Rm(t,-1)和V(t,0)的值，这两个值也没办法求呀。

另外，原贴Rm似乎还有导数项。

不过，这种问题虽然繁琐，但一点一点解应该还是没有问题的吧。

wanglu

看了一下，问题描述的不是很清楚，当k=0是时，v(t,k)为v(t,0)，也就是要求第一、二两个式子Rm(t,-1)和V(t,0)的值，这两个值也没办法求呀。

另外，原贴Rm似乎还有导数项。

不过，这种问题虽然繁琐，但一点一点解 ...
messenger 发表于 2009-12-26 20:25

一时疏忽，问题图片贴错了，已经修改了。另外，V(t,m)对t求一阶导数。
matlab求导比较方便，Forcal中还没有求导函数，我用[f(x+h)-f(x)]/h代替了。

linhengcan126

这个我也在写啊，哈哈，但是很难啊

messenger

大致算了一下，这个题目如果迭代算在技术上没什么难度，只是比较麻烦，只要仔细一点就行了。

Matlab的下标不能从零开始，所以比较麻烦，因为下标不能为零，所以下标要都加1，有的还要减1。

下面是代码，除了下标可能不对（懒得检查了，反正又不是我的问题），其他的技术问题都是对的。

算了一下t=10，m=10的情况，只是几秒钟。

1. 
   
2. 
   
3. syms t h;
   
4. v=sym(zeros(1,11));
   
5. rm=sym(zeros(1,11));
   
6. h=1;
   
7. v(1)=t;
   
8. v(2)=1/3*h*t^3;
   
9. m=10;
   
10. for k1=2:m;
    
11. k=k1+1;
    
12. rm1=diff(v(k-1),t);
    
13. rm2=sym(0);
    
14. 
    
15. for j=1:k1
    
16. temp=v(j)*v(k-j)-(1-k1+1);
    
17. rm2=rm2+temp;
    
18. end
    
19. 
    
20. rm(k-1)=rm1+rm2;
    
21. 
    
22. v(k)=v(k-1)+h*int(rm(k-1),t);
    
23. end
    
24. 
    
25. main=sum(v);
    
26. t=10;
    
27. subs(main,t)
    

复制代码

linhengcan126

thanks 哈哈！

wanglu

看了messenger 版主的代码，虽然不懂matlab的代码，但有些明白了。这个问题用数值解法没有优势，完全可以由解析法解决，而且有很高的效率。

我手算试了一下，用递推公式可以得到v(t,m)和Rm(t,m)的公式，微分式和积分式也可得到。但手算工作量太大，也容易出错。这正好发挥matlab（或mmtc、maple）的优势，先获得解析解，然后计算，就比较简单了。

不知我的理解对不对？

我将这段代码发给原贴主看一下。替他先谢谢messenger 版主了。

另外，另一个问题，放在maple板块：一个含积分的等式优化题
位置：http://forum.simwe.com/thread-910945-1-1.html
是否有解析解呢？

wanglu

不过，关于这个问题我有个预感，是不是当m大到一定数时matlab就算不出了？
这个似乎不是时间问题。
因为随m的增大，符号解的表达式将很长，matlab似乎不能编译太长的表达式?
大家验证一下？

messenger

嗯，是的，这个问题有些碰巧的成份。因为迭代公式都是多项式，积分、求导都比较方便，所以正好可以发挥解析解法的优势。

不过，觉得这个问题如果用数值解法，用数值方法求积分、求导，应该也可以，不过要比解析解法麻烦一些。最重要的是解析解法对自变量 t 没有要求，而数值解法要给出具体的 t ，这也会增加数值积分和求导的难度。

看了messenger 版主的代码，虽然不懂matlab的代码，但有些明白了。这个问题用数值解法没有优势，完全可以由解析法解决，而且有很高的效率。

我手算试了一下，用递推公式可以得到v(t,m)和Rm(t,m)的公式，微分式和积分式也可得到。但手算工作量太大，也容易出错。这正好发挥matlab（或mmtc、maple）的优势，先获得解析解，然后计算，就比较简单了。

不知我的理解对不对？
wanglu 发表于 2009-12-27 17:44

messenger

如果表达式太长的话，可能运算时间较长。不过，表达式的长度限制倒不是问题，因为此问题是加法运算，可以设中间变量，将一个较长的表达式拆分。比如一个100项相加的表达式，代换成2个变量相加，每一个变量为50项相加。

不过，关于这个问题我有个预感，是不是当m大到一定数时matlab就算不出了？
这个似乎不是时间问题。
因为随m的增大，符号解的表达式将很长，matlab似乎不能编译太长的表达式?
大家验证一下？
wanglu 发表于 2009-12-27 18:14

wanglu

如果表达式太长的话，可能运算时间较长。不过，表达式的长度限制倒不是问题，因为此问题是加法运算，可以设中间变量，将一个较长的表达式拆分。比如一个100项相加的表达式，代换成2个变量相加，每一个变量为50项相加 ...
messenger 发表于 2009-12-27 19:29

明白了。

另外，积分式中那个指数2k中的k取0是一个特殊情况，原贴主为了简化计算，将k取0了。若k不为0，而是3.2或者是其他数值，解析解法还能求解吗？

messenger

k不为零应该也没问题，原来被积函数Rm为多项式，再乘一个tau^k，还是多项式，这对积分来说没有什么新困难。

另外，积分式中那个指数2k中的k取0是一个特殊情况，原贴主为了简化计算，将k取0了。若k不为0，而是3.2或者是其他数值，解析解法还能求解吗？
wanglu 发表于 2009-12-27 20:01

linhengcan126

这个一个递归或迭代求解问题目的用来干什么的呢？？

linhengcan126

h=1;
v(1)=t;
v(2)=1/3*h*t^3;
其中的作用是干什么用的呢？？

messenger

v(1)、v(2)是初值呀，递归或递推都要有初值呀。比如x(i+2)=x(i+1)+x(i)，要知道x(1)和x(2)的值才能知道x(3)的值，进而推得其他x的值。

h=1;
v(1)=t;
v(2)=1/3*h*t^3;
其中的作用是干什么用的呢？？
linhengcan126 发表于 2009-12-27 21:01

wanglu

k不为零应该也没问题，原来被积函数Rm为多项式，再乘一个tau^k，还是多项式，这对积分来说没有什么新困难。


messenger 发表于 2009-12-27 20:50

谢谢了！
看来解析解法和数值解法各有所长，需要用其所长啊。

wanglu

h=1;
v(1)=t;
v(2)=1/3*h*t^3;
其中的作用是干什么用的呢？？
linhengcan126 发表于 2009-12-27 21:01

这个我也不清楚，请查看原帖。
原帖位置：http://www.ilovematlab.cn/thread-60022-1-1.html

linhengcan126

有时候我在怀疑，我们学习这么深，这些函数的真正意义是干什么的呢？？

messenger

递推和递归在实际问题中还是很常见的

有时候我在怀疑，我们学习这么深，这些函数的真正意义是干什么的呢？？
linhengcan126 发表于 2009-12-27 21:22

linhengcan126

能不能举个例子，我感觉不到啊，thanks!