龙格库塔算法

全部家当了，有了还可以追加，求大神拉一把 TT

1. ##runge-kutta
   
2. 
   
3. Runge_Kutta <- function(f,b,c,A,y0,x0,h,N){
   
4.     stopifnot( length(c)== nrow(A))
   
5.     x=x0
   
6.     y=y0
   
7.     message(sprintf( 'initial condition, x=%f, y=%f',x,y))
   
8.     for(n in 1:N){
   
9.                                         #cal k1
   
10.         K=h*f(x,y)
    
11.                                         #cal k2...kN
    
12.         for(s in 2:length(c)){
    
13.             k= h*f(x+c[s]*h, y+ A[s, 1:(s-1) ]*K )
    
14.             K=c(K,k)
    
15.         }
    
16. 
    
17.                                         #update x,y
    
18.         y=y+ sum(b*k)
    
19.         x=x+ h
    
20.         message(sprintf( 'iteration %d, x=%f, y=%f',n,x,y))
    
21. 
    
22.     }
    
23.     yN=y
    
24. 
    
25.     return(yN)
    
26. }
    
27. 
    
28. 
    
29. b= c(1/6,1/3,1/3,1/6)
    
30. c= c(0, 1/2,1/2,1)
    
31. A= matrix(c(0,0,0,0,1/2,0,0,0,0,1/2,0,0,0,0,1,0), nrow=4,ncol=4,byrow=T)
    
32. f=function(x,y){-2*y*x^2}
    
33. 
    
34. yN=Runge_Kutta(f,b,c,A,y0=1,x0=0,h=.001,N=1000)
    
35. 
    
36. yE=exp(-2/3* 1^3)
    
37. message(sprintf('Runge_Kutta method result at  x=1 is %f; exactly result is %f.',yN,yE))

复制代码

发重了，删了吧~

马甲1号 发表于 2016-5-10 13:43
发重了，删了吧~

谢谢大神[cry]