Swift 3中的Lotka-Volterra模型

Question

我正在尝试实现一个rk4函数来求解2个微分方程。我有这样的代码，实现了龙格库塔4方法：

//RK4 method 
    func rk4_func(y_array: [Double], f_array: [(([Double], Double) -> Double)], t_val: Double, h_val: Double) -> [Double] { 

     let length = y_array.count 

     let t_half_step = t_val + h_val/2.0 
     let t_step = t_val + h_val 

     var k1 = [Double](repeating: 0.0, count: length) 
     var k2 = [Double](repeating: 0.0, count: length) 
     var k3 = [Double](repeating: 0.0, count: length) 
     var k4 = [Double](repeating: 0.0, count: length) 
     var w = [Double](repeating: 0.0, count: length) 
     var result = [Double](repeating: 0.0, count: length) 

     for i in 0...length { 
      k1[i] = h_val * f_array[i](y_array, t_val) 
      w[i] = y_array[i] + k1[i]/2.0 
     } 

     for i in 0...length { 
      k2[i] = h_val * f_array[i](w, t_half_step) 
      w[i] = y_array[i] + k2[i]/2.0 
     } 

     for i in 0...length { 
      k3[i] = h_val * f_array[i](w, t_half_step) 
      w[i] = y_array[i] + k3[i] 
     } 

     for i in 0...length { 
      k4[i] = h_val * f_array[i](w, t_step) 
     } 

     for i in 0...length { 
      result[i] = y_array[i] + (k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i])/6.0 
     } 

     print(result) 
     return result; 

    } 

但现在我需要真正使用它，这是我感到困惑的部分。如果有人对微分方程有数值计算方法的经验，那将会有所帮助。

1. 需要哪些数组来提供此功能？

2. t_val参数代表什么？这是最大时间吗？

3. 输出如何“解决”方程？

4. 输出给我什么？

5. 在行k1[i] = h_val * f_array[i](y_array, t_val)，f_array[i](y_array, t_val)是什么意思？是否说对于f_array的i-th值，找到y_array对应的i-th值？那么t_val是什么意思？

仅供参考，以下是需要解决的2个微分方程。上下文是我试图数值求解这些Lotka-Volterra模型来绘制Xcode（Swift 3.x）中的时间序列和相位空间图。

Answer 1

1. y是当前状态（如双阵列实现）的向量。 f_array是函数指针doty = f_array(y,t)。

2. t_val是当前状态的时间，h_val是时间步长。的rk4_func

3. 一个呼叫执行从t_val到t_val+h_val和

4. 返回新的状态，y_next = rk4_func(y, f_array, t, h)的时间步长。

5. 人们将不​​得不学习语言内部。希望，也就是说，要使代码正常工作，第一次调用f_array[0](y_array, t_val)将计算完整的向量/数组值，并进一步调用提取缓存结果的组件。

---

截至https://github.com/pdemarest/swift-rk4发现原始代码是其实现RK4和语言标准外的日期严重缺乏。截至https://swift.sandbox.bluemix.net/测试工作版本是

import Foundation 


func RK4step(y: [Double], f: ([Double], Double) -> [Double], t: Double, h: Double) -> [Double] { 

    let length = y.count 

    var w = [Double](repeating: 0.0, count: length) 
    var result = [Double](repeating: 0.0, count: length) 

    let k1 = f(y,t) 
    assert(k1.count == y.count, "States and Derivatives must be the same length") 
    for i in 0..<length { w[i] = y[i] + 0.5*h*k1[i] } 
    let k2 = f(w, t+0.5*h) 
    for i in 0..<length { w[i] = y[i] + 0.5*h*k2[i] } 
    let k3 = f(w,t+0.5*h) 
    for i in 0..<length { w[i] = y[i] + h*k3[i] 
    } 
    let k4 = f(w,t+h) 

    for i in 0..<length { 
     result[i] = y[i] + (k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i])*h/6.0 
    } 

    return result; 
} 



func test_exp(){ 
    // Integrate: y' = y 
    // y_0 = 1.0 
    // from 0 to 2.0 

    var y = [1.0] 

    func deriv (y: [Double], t: Double) -> [Double] { 
     return [ y[0] ] 
    } 

    var t = 0.0 
    let h = 0.1 


    while t < 2.0 { 
     y = RK4step(y:y, f:deriv, t:t, h:h) 
     t += h 
     print("t: \(t), y: \(y[0]) exact: \(exp(t))\n") 
    } 

    let exact = exp(2.0) 
    let error = abs(y[0] - exact) 

    assert(error < pow(h, 4.0)) 
    print("y: \(y[0]) exact: \(exact) error: \(error)\n") 

} 


print("testing...\n") 

test_exp() 

对于沃尔泰拉 - 洛特卡动力学一个将不得不改变

var y = [150.0, 5.0] 

    let a = 5.0 
    let b = 1.0 
    let eps = 0.1 
    let m = 5.0 

    func deriv (y: [Double], t: Double) -> [Double] { 
     let p = y[0] 
     let q = y[1] 
     return [ a*p-b*p*q, eps*b*p*q - m*q ] 
    } 

与妥善固定全局常量a,b,eps,m和二维初始值。在需要时添加打印语句。