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

linhengcan126

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


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

如果k不等于0 那如何修改上面的程序呢？？thanks!

messenger

这与专业有关系，有的专业用的多，尤其是循环现象较多的专业，比如振动问题的滞洄现象，Simulink仿真中也常用到，只要是下一个状态依赖上一个状态，就会有递推或递归。

一个比较简单的例子，就是经济增长，每年的增长比上一年增长10%，1980年是1000亿元，问2010年是多少元？当然这一个问题有公式可算，但如果复杂一点的，就没有公式来算了，就要用递推或递归来算了，而且那个公式的推导也是递推出来的。一般递推比递归用的多，毕竟递归不太符合正常思维。

能不能举个例子，我感觉不到啊，thanks!
linhengcan126 发表于 2009-12-27 21:42

messenger

将v(k)=v(k-1)+h*int(rm(k-1),t);改为v(k)=v(k-1)+h*int(t^(k-1)*rm(k-1),t);

如果k不等于0 那如何修改上面的程序呢？？thanks!
linhengcan126 发表于 2009-12-27 22:02

linhengcan126

thanks messenger ！哈哈！

linhengcan126

这与专业有关系，有的专业用的多，尤其是循环现象较多的专业，比如振动问题的滞洄现象，Simulink仿真中也常用到，只要是下一个状态依赖上一个状态，就会有递推或递归。

一个比较简单的例子，就是经济增长，每年的 ...
messenger 发表于 2009-12-27 22:14

一个比较简单的例子，就是经济增长，每年的增长比上一年增长10%，1980年是1000亿元，问2010年是多少元？当然这一个问题有公式可算，但如果复杂一点的，就没有公式来算了，就要用递推或递归来算了，而且那个公式的推导也是递推出来的

如果像这个题目一样，那么多函数，如何才能推导和递推一个那么复杂的公式呢？？

linhengcan126

所以我就不明白，为什么会有那么复杂的公式出现，是用来模拟什么事物的呢？？有那么复杂吗？

messenger

,
等你遇到了自然就知道了，其实我遇到的复杂的递归/递推也不多，就本版问题来看，涉及这一问题的也不多。

其实，这种问题虽然复杂，但解决起来也不困难，只要仔细一些，用循环都能解决。

lengyunfeng

所以我就不明白，为什么会有那么复杂的公式出现，是用来模拟什么事物的呢？？有那么复杂吗？
linhengcan126 发表于 2009-12-27 22:43

现实生活中的问题、信息不是那么纯粹、直接的，要用数学式子来对这些信息进行重现已经是很难得的了，更不用说想只用一个式子表达清楚一个问题。数学发展到今天，虽然已经解决了部分非确定性问题，但还有很多前沿科学问题是不可能完全解，比如灰色理论、模糊数学。我们最大的期望是能用一个个的数学式子表达清楚每一句话的意思，然后再图将其化简的办法。你所遇到的不应该只是初高中的应用题吧？朋友，是时候该考虑一下向银行借贷买房，如果银行每天贷款利率都在变，你一个月还几千块钱，那得啥时候才能还完哦～～～

xl0418

28# lengyunfeng

其实这个迭代方程式一个称作Homotopy Analysis Method的求解非线性微分方程的算法的一部分步骤。我的目的是想验证这个方法的精确性，以及应用到另一个非线性微分方程求解上。

wanglu

原贴主到了，这样讨论问题比较方便。
问题图片我又更新了，积分式中的τ实际上应是t？

我优化了Forcal的数值算法，主要是为了提高速度，方法：保存中间变量，用插值法近似取值。

但效果不算太好。求解精度取1e-3，算t＝10，m＝5需要几分钟的时间，其他的没有尝试。

1. //SetTalk[false];
   
2. !SetRealStackMax(1000);
   
3. Rm(t,m:j,s,a,b,h:hVm,iRm,tmax)=
   
4. {
   
5.   iRm++, h=0.0005, if[t>=h, h=-0.0005],
   
6.   sys::call(2,hVm :t+h,m : &a),  //计算s=Vm(t+h,m)
   
7.   sys::call(2,hVm :t,m : &b),    //计算s=Vm(t,m)
   
8.   s=[a-b]/h,
   
9.   j=0,(j<=m).while{
   
10.     sys::call(2,hVm :t,j : &a), sys::call(2,hVm :t,m-j : &b),  //计算a=Vm(t,j),   b=Vm(t,m-j)
    
11.     s=s+a*b,
    
12.     j++
    
13.   },
    
14.   s-[1-which(m<1,0,1)] //最后返回这个值。看这句有没有问题：取m<1，而不是m<=1
    
15. };
    
16. //Fkappa将返回s的值。mm是模块变量，因调用Rm(t,mm-1)时mm将改变，故用a预先保存，函数返回前恢复mm的值
    
17. Fkappa(t:a,s,i,kappa,x1,x2,y1,y2:mm,_t,_m,tmax,VmArray,iFkappa)=
    
18. { iFkappa++, kappa=0, a=mm,
    
19.   i=t/_t*tmax,
    
20.   if{VmArray.XSLSF::getrai[a+_m+1,i,0] & VmArray.XSLSF::getrai[a+_m+1,i+1,0],      //获得Fkappa(t,m)
    
21.      x1=VmArray.XSLSF::getrai(a+_m+1,i,1), y1=VmArray.XSLSF::getrai(a+_m+1,i,2),
    
22.      x2=VmArray.XSLSF::getrai(a+_m+1,i+1,1), y2=VmArray.XSLSF::getrai(a+_m+1,i+1,2),
    
23.      return[y1+(y2-y1)*(t-x1)/(x2-x1)]
    
24.   },
    
25.   s=t^(2*kappa)*Rm(t,mm-1), mm=a,
    
26.   i=t/_t*tmax, VmArray.XSLSF::setrai[a+_m+1,i,0,1], VmArray.XSLSF::setrai[a+_m+1,i,1,t], VmArray.XSLSF::setrai[a+_m+1,i,2,s], //保存Fkappa(t,m)
    
27.   s
    
28. };
    
29. Vm(t,m:h,kappa,s,i,x1,x2,y1,y2:VmArray,_t,tmax,mm,hFkappa,iVm)=
    
30. {
    
31.   iVm++, h=-1, kappa=0,
    
32.   if[m<=0,return(t)], if[m==1,return(1/3*h*t^3)],
    
33.   i=t/_t*tmax,
    
34.   if{VmArray.XSLSF::getrai[m,i,0] & VmArray.XSLSF::getrai[m,i+1,0],      //获得Vm(t,m)
    
35.      x1=VmArray.XSLSF::getrai(m,i,1), y1=VmArray.XSLSF::getrai(m,i,2),
    
36.      x2=VmArray.XSLSF::getrai(m,i+1,1), y2=VmArray.XSLSF::getrai(m,i+1,2),
    
37.      return[y1+(y2-y1)*(t-x1)/(x2-x1)]
    
38.   },
    
39.   s=which(m<=1,0,1)*Vm(t,m-1)+h*{mm=m, XSLSF::fpqg[hFkappa,0,t,1e-3]},         //XSLSF::fpqg即徐士良算法中的连分式积分法
    
40.   i=t/_t*tmax, VmArray.XSLSF::setrai[m,i,0,1], VmArray.XSLSF::setrai[m,i,1,t], VmArray.XSLSF::setrai[m,i,2,s], //保存Vm(t,m)
    
41.   s
    
42. };
    
43. pert(t)=t-1/3*t^3+2/15*t^5-17/315*t^7;
    
44. main(t,m:k,s:VmArray,_t,_m,tmax,hFkappa,hVm,iRm,iFkappa,iVm)=
    
45. {
    
46.   iRm=0, iFkappa=0, iVm=0,
    
47.   _t=t, _m=m, tmax=20000,  //tmax>=10000，越大越精确
    
48.   VmArray=new[rtoi(real_s),rtoi(m+m+2),rtoi(tmax+2),rtoi(3)],
    
49.   k=0,(k<=m+m+1).while{
    
50.     s=0,(s<=tmax).while{
    
51.       VmArray.XSLSF::setrai[k,s,0,0],
    
52.       s++
    
53.     },
    
54.     k++
    
55.   },
    
56.   hFkappa=HFor("Fkappa"), hVm=HFor("Vm"),  //预先获得函数Fkappa和Vm的句柄，函数调用时将用到
    
57.   s=0,k=0,(k<=m).while{
    
58.     s=s+Vm(t,k),
    
59.     k++
    
60.   },
    
61.   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},
    
62.   delete[VmArray],
    
63.   s  //返回s
    
64. };
    
65. main[10,3];
    
66. main[10,5];

复制代码

m＝3时的验证公式：pert(t)=t-1/3*t^3+2/15*t^5-17/315*t^7;

从t＝10，m＝3的验证结果看，应该是准确的（t＝10，m＝5时没有验证公式）：

结果＝-526669.7208717888       , 验证pert(t)＝-526672.5396825398       , 函数Rm调用次数74669   , 函数Fkappa调用次数149094  , 函数Vm调用次数451097  
结果＝-1206839696.504438       , 验证pert(t)＝-526672.5396825398       , 函数Rm调用次数334056  , 函数Fkappa调用次数68015611, 函数Vm调用次数2799112

messenger

觉得第2、3个公式中有关Rm和V的公式里，t 和 m 应该是 tau 和 k。 t 和 m 都是外层变量。

另外，如果如lz所说，这个问题是由解其他方程过程中转换过来的，最好还是用数值解，这个运行效率可能好一些，不用做符号与数值的转换。

wanglu

原贴主的原图我看起来也比较费劲，然后我改写了一下。目前图的描述与Forcal代码完全一致，且结果与原贴主给出的验证公式一致。是不是还有问题，需要等原帖主更正了。

刚刚又算了一下t＝10，m＝6，用了10多分钟，函数调用约2亿次了，不能再试了。

结果＝-11972514105864.28       , 验证pert(t)＝-526672.5396825398       , 函数Rm调用次数838630  , 函数Fkappa调用次数190977727, 函数Vm调用次数7001144
-11972514105864.28