0

J'essaye d'implémenter une fonction rk4 pour résoudre 2 équations différentielles. J'ai ce code qui implémente la méthode Runge Kutta 4:Lotka-Volterra Modèles dans Swift 3

//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; 

    } 

Mais maintenant, je dois effectivement utilisation, ce qui est la partie que je suis confus au sujet. Si quelqu'un a de l'expérience en calcul numérique des solutions aux équations différentielles, cela aiderait.

  1. De quels tableaux ai-je besoin pour alimenter cette fonction?

  2. Que représente l'argument t_val? Est-ce un temps maximum?

  3. Comment la sortie "résout" l'équation?

  4. Que me donne la sortie? Dans la ligne k1[i] = h_val * f_array[i](y_array, t_val), que signifie f_array[i](y_array, t_val)? Est-ce que cela signifie que pour la valeur i-th de f_array, trouvez la valeur i-th correspondante pour y_array? Alors qu'est-ce que le t_val signifie là?

Pour référence, voici les 2 équations différentielles à résoudre. Le contexte est que j'essaie de résoudre numériquement ces modèles de Lotka-Volterra pour tracer une série temporelle et un tracé d'espace de phase dans Xcode (Swift 3.x).

enter image description here

enter image description here

+0

Dans l'appel 'k2 [i] = h_val * f_array [i] (w, t_half_step)', le résultat est-il calculé comme un vecteur lors de la première exécution de la boucle? Si ce n'est pas le cas, la valeur changée de 'w [i]' ferait des ravages avec la méthode. – LutzL

+0

Je pense que vous voulez dire 'array', pas vectoriel. Swift utilise 'arrays'. Et oui, toutes les variables 'k' sont' arrays'. – loltospoon

+0

La question est de savoir à quel moment 'f_array [1]' est évalué? Pour que tout se passe bien, il devrait être simultanément avec 'f_array [0]'. Et cela rendrait le code plus simple si cela pouvait être fait avec des vecteurs doubles (impliquant des opérations vectorielles) au lieu de doubles tableaux nus. – LutzL

Répondre

0
  1. y est le vecteur de l'état actuel (mis en œuvre en tant que matrice double). f_array est un pointeur de fonction vers une fonction doty = f_array(y,t).

  2. t_val est le temps pour l'état actuel, h_val est le pas de temps.

  3. Un appel rk4_func effectue le pas de temps de t_val à t_val+h_val et

  4. retourne le nouvel état, y_next = rk4_func(y, f_array, t, h).

  5. Il faudrait étudier les composants internes de la langue. Heureusement, pour que le code fonctionne correctement, le premier appel de f_array[0](y_array, t_val) calcule le résultat vectoriel/tableau-complet et d'autres appels extraient simplement les composants du résultat mis en cache.


Le code original que l'on trouve à https://github.com/pdemarest/swift-rk4 est gravement déficiente dans sa réalisation RK4 et hors jour dans les normes linguistiques.Une version de travail testé à https://swift.sandbox.bluemix.net/ est

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() 

Pour la dynamique Volterra-Lotka on aurait pu changer

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 ] 
    } 

avec des constantes globales correctement fixes a,b,eps,m et une valeur initiale en deux dimensions. Ajoutez des instructions d'impression si nécessaire.

+0

En ce qui concerne # 1, pourquoi ai-je besoin de passer dans 'f_array: [(([Double], Double) -> Double)]' au lieu de simplement 'f_array: [Double]'? – loltospoon

+0

Parce que c'est le pointeur/référence au sous-programme qui implémente la fonction de l'ODE 'y '= f (y, t)', pas un tableau constant. – LutzL

+0

Ok. Donc, dans 'f_array [i] (y_array, t_val)' cela signifie-t-il que nous passons '(y_array, t_val)' en arguments? – loltospoon