0

I'm trying to get more acquainted with Clojre so I decided to do my Runge Kutta integrator project in it. However I'm having problems working with the immutable nature of the let statement. I want to evaluate 8 variables in each iteration of the loop and use them to recurse through it until my loop is finished.

The way I understand it, since my recur is inside the let's scope, my k's and l's won't be overwritten with each recursion. I'm looking for a more idiomatic way to recurse through my integrator.

(loop [psin 0  Xin 1 xn dx] ;initial values
  (if (> xn 1)
    psin
    (let [k1 (evalg  xn psin) ;define 4 variables to calculate next step of Xi, psi
          l1 (evalf  xn Xin)  ;evalf and evalg evaluate the functions as per Runge Kutta
          k2 (evalg (+ (* 0.5 dx) xn) (+ (* 0.5 l1) psin))
          l2 (evalf (+ (* 0.5 dx) xn) (+ (* 0.5 k1) Xin))
          k3 (evalg (+ (* 0.5 dx) xn) (+ (* 0.5 l2) psin))
          l3 (evalf (+ (* 0.5 dx) xn) (+ (* 0.5 k2) Xin))
          k4 (evalg (+ dx xn) (+ l3 psin))
          l4 (evalf (+ dx xn) (+ k3 Xin))]
      (do
        (let [Xinew (+ Xin (* (/ dx 6) (+ k1 k4 (* 2 k3) (* 2 k2))) )
              psinew (+ psin (* (/ dx 6) (+ l1 l4 (* 2 l2) (* 2 l3) )))]
          (println k1)
          (recur psinew Xinew (+ dx xn)))))))

Many thanks! Looking forward to getting more acquainted with clojure:)

apmechev
  • 75
  • 1
  • 8
  • I don't think your understanding is correct - the `let` bindings will be re-evaluated on each loop iteration. Did you try it? – Alex Oct 29 '14 at 15:17

1 Answers1

0

First, as @Alex has commented, the recur does evaluate your ks and ls afresh on every iteration of the loop. A recur refers to its immediately enclosing function or loop form, in this case the latter.


The (println k1) suggests that what you are looking for is a sequence of k1s as you iterate round the loop.

Clojure has the concept of a lazy sequence: a potentially endless sequence of values that are

  • calculated when called for
  • forgotten if/when no longer accessible (garbage collected).

You should

  • frame the solution as a lazy-sequence of k1s and - while we're at it -
  • make your evalg and evalf functions arguments to your runge-kutta function.

We could build our lazy sequence from scratch, using explicit recursion wrapped in the lazy-seq macro, but there is a whole library of sequence functions to hand, one of which will usually express your intent. These can be faster than what you or I would write, but will still be much slower than your loop.

A handy function here is iterate. We can systematically convert your loop to use it as follows:

  • Convert the several loop arguments into a single destructured vector argument to a local function (which I've called step).
  • Include any result you want, in this case k1 in the argument list.
  • Rewrite the iteration to return what were the next set of arguments to the loop as a vector.
  • Call map on the resulting sequence to get the data you want.

My best guess is the following:

(defn runge-kutta [evalg evalf dx]
  (letfn [(step [[k1 psin  Xin xn]]
                (let [k1 (evalg  xn psin)
                      l1 (evalf  xn Xin)
                      k2 (evalg (+ (* 0.5 dx) xn) (+ (* 0.5 l1) psin))
                      l2 (evalf (+ (* 0.5 dx) xn) (+ (* 0.5 k1) Xin))
                      k3 (evalg (+ (* 0.5 dx) xn) (+ (* 0.5 l2) psin))
                      l3 (evalf (+ (* 0.5 dx) xn) (+ (* 0.5 k2) Xin))
                      k4 (evalg (+ dx xn) (+ l3 psin))
                      l4 (evalf (+ dx xn) (+ k3 Xin))
                      Xinew (+ Xin (* (/ dx 6) (+ k1 k4 (* 2 k3) (* 2 k2))) )
                      psinew (+ psin (* (/ dx 6) (+ l1 l4 (* 2 l2) (* 2 l3) )))]
                  [k1 psinew Xinew (+ dx xn)]))]
    (map first (iterate step [nil 0 1 dx]))))

There may well be errors in this, as I'm groping in the dark.

The sequence is endless. You can stop it

  • by returning the xns too, and wrapping the result in a take-while or
  • working out how may iterations you want - it looks like (long (/ dx)), though there may be a one-off error there. Then just use nth or take to get what you want.

Let us know how you get on.

Thumbnail
  • 13,293
  • 2
  • 29
  • 37