【独家发布】R语言计算费根鲍姆常数

R语言计算费根鲍姆常数
完成标准：
1、完全通过R语言实现，不调用其他语言外部代码；
2、算法部分请注释；
3、计算到小数点后30位。

1. library(Rmpfr)
   
2. 
   
3. # 具体参考 [url]http://keithbriggs.info/documents/how-to-calc.pdf[/url]
   
4. 
   
5. maxIt = 13  
   
6. maxItJ = 10
   
7. prec=30*log2(10) #小数点后30位
   
8. a_1 = Rmpfr::mpfr(1.0, prec)
   
9. a_2 = Rmpfr::mpfr(0.0, prec)
   
10. d_1 = Rmpfr::mpfr(3.2, prec)
    
11. 
    
12. res = NULL
    
13. 
    
14. for (i in 2:maxIt) {
    
15.   a = a_1 + (a_1 - a_2)/d_1
    
16.   for (j in 1:maxItJ) {
    
17.     x = 0
    
18.     y = 0
    
19.     for (k in 1:2^i) {
    
20.       y = 1 - 2*y*x
    
21.       x = a - x^2
    
22.     }
    
23.     a = a - x/y
    
24.   }
    
25.   d = (a_1 - a_2)/(a-a_1)
    
26.   #print(c(i,a,d))
    
27.   d_1 = d
    
28.   a_2 = a_1
    
29.   a_1 = a
    
30.   res = c(res,d)
    
31. }
    
32. 
    
33. print(res)
    
34. 
    
35. # 'mpfr1' 4.669201545780906707452765346774

复制代码

zhou1_20 发表于 2021-6-23 15:15

谢谢zhou1_20老师。
您提供的代码已经是我找到的所有代码中最接近的了。所以，虽然与设想有异，我还是准备结贴了。

以下只是请教您：
不知道那个一直在算的结果（比如我附的截图中的官方结果），采用的是什么方式。
另外的问题
前一段时间，尝试了verhulst模型的分叉图，感觉只有第1、2、3次分叉，最多第4次分叉，用肉眼能分辨出来，而用模型的离散方程做模拟，增长率r大于3.1时，曲线直接就崩溃了。因此我猜这个r只能是不断接近某个极限，只是不知道这个极限是什么决定的。

我再啃一下您提到的那篇文献。
明天有时间，我再试试您的代码，看看增加迭代次数会是什么结果。

再次感谢。

提示：费根鲍姆常数的绘图和求解是两种编程思路

可参考百度文库的这一篇以及我上面提到的pdf
https://wenku.baidu.com/view/d5d ... PrWF5MeE%2FaA%3D%3D

1. logistic.map <- function(r, x, N, M){
   
2.   ## r: bifurcation parameter
   
3.   ## x: initial value
   
4.   ## N: number of iteration
   
5.   ## M: number of iteration points to be returned
   
6.   z <- 1:N
   
7.   z[1] <- x
   
8.   for(i in c(1:(N-1))){
   
9.     z[i+1] <- r *z[i]  * (1 - z[i])
   
10.   }
    
11.   ## Return the last M iterations
    
12.   z[c((N-M):N)]
    
13. }
    
14. 
    
15. ## Set scanning range for bifurcation parameter r
    
16. my.r <- seq(2.5, 4, by=0.003)
    
17. system.time(Orbit <- sapply(my.r, logistic.map,  x=0.1, N=1000, M=300))
    
18. 
    
19. Orbit <- as.vector(Orbit)
    
20. 
    
21. r <- sort(rep(my.r, 301))
    
22. 
    
23. plot(Orbit ~ r, pch=".")

复制代码

zhou1_20 发表于 2021-6-24 09:14
提示：费根鲍姆常数的绘图和求解是两种编程思路

可参考百度文库的这一篇以及我上面提到的pdf

感谢！非常感谢！