1. import os
  2. import strutils
  4. let params = commandLineParams()
  6. let n = parseInt(params[1])
  7. let h = parseFloat(params[2])
  8. let steps = parseInt(params[3])
  10. var x0 = parseFloat(params[4])
  11. var y0 = parseFloat(params[5])
  12. var z0 = parseFloat(params[6])
  14. let sigma = parseFloat(params[7])
  15. let rho = parseFloat(params[8])
  16. let beta = parseFloat(params[9]) / parseFloat(params[10])
  18. var x = newSeq[float64](n + 1)
  19. var y = newSeq[float64](n + 1)
  20. var z = newSeq[float64](n + 1)
  22. proc t_horner (j: openArray[float64], h: float64): float64 =
  23. for i in countdown(j.len - 1, 0):
  24. result = result * h + j[i]
  26. proc t_prod (u, v: openArray[float64], k: int): float64 =
  27. for j in 0..k:
  28. result += u[j] * v[k - j]
  30. echo(x0, ' ', y0, ' ', z0, ' ', 0.0)
  31. for step in 1..steps+1:
  32. x[0] = x0
  33. y[0] = y0
  34. z[0] = z0
  35. for k in 0..n-1:
  36. x[k + 1] = sigma * (y[k] - x[k]) / float64(k + 1)
  37. y[k + 1] = (rho * x[k] - y[k] - t_prod(x, z, k)) / float64(k + 1)
  38. z[k + 1] = (t_prod(x, y, k) - beta * z[k]) / float64(k + 1)
  39. x0 = t_horner(x, h)
  40. y0 = t_horner(y, h)
  41. z0 = t_horner(z, h)
  42. echo(x0, ' ', y0, ' ', z0, ' ', float64(step) * h)

Lorenz ODE Solver

nim compile --run test.nim 16 10 .01 10 -15.8 -17.48 35.64 10 28 8 3