runge kutta第二种方法的python代码是什么？

这个代码是否可用？runge kutta第二种方法的python代码是什么？

def rKN(x, fx, n, hs): 
    k1 = [] 
    k2 = [] 
    k3 = [] 
    k4 = [] 
    xk = [] 
    for i in range(n): 
     k1.append(fx[i](x)*hs) 
    for i in range(n): 
     xk.append(x[i] + k1[i]*0.5) 
    for i in range(n): 
     k2.append(fx[i](xk)*hs) 
    for i in range(n): 
     xk[i] = x[i] + k2[i]*0.5 
    for i in range(n): 
     k3.append(fx[i](xk)*hs) 
    for i in range(n): 
     xk[i] = x[i] + k3[i] 
    for i in range(n): 
     k4.append(fx[i](xk)*hs) 
    for i in range(n): 
     x[i] = x[i] + (k1[i] + 2*(k2[i] + k3[i]) + k4[i])/6 
    return x

来源

2016-02-07 Amin Mithun

考虑是否行得通？你怎么知道的？ –

查看https://stackoverflow.com/questions/35071393/runge-kutta-code-not-converging-with-builtin-method（或侧栏上的任何“相关”链接）以获得更紧凑的RK4代码。 – LutzL

与numpy似乎更具可读性：

代码从http://www.math-cs.gordon.edu/courses/mat342/python/diffeq.py

def rk2a(f, x0, t): 
    """Second-order Runge-Kutta method to solve x' = f(x,t) with x(t[0]) = x0. 

    USAGE: 
     x = rk2a(f, x0, t) 

    INPUT: 
     f  - function of x and t equal to dx/dt. x may be multivalued, 
       in which case it should a list or a NumPy array. In this 
       case f must return a NumPy array with the same dimension 
       as x. 
     x0 - the initial condition(s). Specifies the value of x when 
       t = t[0]. Can be either a scalar or a list or NumPy array 
       if a system of equations is being solved. 
     t  - list or NumPy array of t values to compute solution at. 
       t[0] is the the initial condition point, and the difference 
       h=t[i+1]-t[i] determines the step size h. 

    OUTPUT: 
     x  - NumPy array containing solution values corresponding to each 
       entry in t array. If a system is being solved, x will be 
       an array of arrays. 

    NOTES: 
     This version is based on the algorithm presented in "Numerical 
     Analysis", 6th Edition, by Burden and Faires, Brooks-Cole, 1997. 
    """ 

    n = len(t) 
    x = numpy.array([ x0 ] * n) 
    for i in xrange(n - 1): 
     h = t[i+1] - t[i] 
     k1 = h * f(x[i], t[i])/2.0 
     x[i+1] = x[i] + h * f(x[i] + k1, t[i] + h/2.0) 

    return x

来源

2016-02-07 20:51:51

runge kutta第二种方法的python代码是什么？

回答

相关问题