I wrote fluid advection code: How to make this more elegant?

I wrote fluid advection code: How to make this more elegant?

[Arne Babenhauserheide]

[-- Attachment #1: Type: text/plain, Size: 4745 bytes --]

Hi,

I just recreated a fluid advection exercise in Guile Scheme and I’m not
quite happy with its readability. Can you help me improve it?

My main gripe is that the math does not look instantly accessible.

The original version was in Python:

    psi[i] - c1*(psi[i+1] - psi[i-1]) + c2*(psi[i+1] - 2.0*psi[i] + psi[i-1])

My port to Scheme looks like this:

    (let ((newvalue (+ (- (psir i)
                          (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
                       (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
                                (psir (- i 1)))))))
      (array-set! psinew newvalue i))


Liebe Grüße,
Arne

Here’s the full code:

#!/bin/sh
# -*- scheme -*-
exec guile -e '(@@ (advection) main)' -s "$0" "$@"
!#

; Copyright (c) 2015 John Burkardt (original Python), 2016 Corinna
; Hoose (adaption) and 2016 Arne Babenhauserheide (pep8 + Scheme
; version).

; License: LGPL, built on the Python version from 2015 John Burkardt
; and Corinna Hoose. License LGPL.

(define-module (advection)
  #:use-module (ice-9 optargs) ; define*
  #:use-module (srfi srfi-1) ; iota
  #:use-module (ice-9 format)
  #:use-module (ice-9 popen))


(define* (fd1d-advection-lax-wendroff #:key (nx 101) (nt 1000) (c 1))
  (let* ((dx (/ 1 (- nx 1)))
         (x (iota nx 0 (/ 1 nx)))
         (dt (/ 1 nt))
         (c1 (* #e0.5 (* c (/ dt dx))))
         (c2 (* 0.5 (expt (* c (/ dt dx)) 2))))
    (format #t "CFL condition: dt (~g) ≤ (~g) dx/c\n" dt (/ dx c))
    (let ((psi (make-array 0 nx))
          (X (make-array 0 nx (+ nt 1)))
          (Y (make-array 0 nx (+ nt 1)))
          (Z (make-array 0 nx (+ nt 1))))
      (let ((psinew (let ((pn (make-array 0 nx)))
                      (let loop ((i 0))
                        (cond ((= i nx)
                               pn)
                              (else
                               (let ((xi (list-ref x i)))
                                 (when (and (<= 0.4 xi) (<= xi 0.6))
                                   (array-set! pn
                                               (* (expt (- (* 10 xi) 4) 2)
                                                  (expt (- 6 (* 10 xi)) 2))
                                               i))
                                 (loop (+ 1 i)))))))))
        (define (psir i) (array-ref psi i))
        (let loop ((j 0))
          (cond
           ((> j nt) #t) ; done
           (else 
            (let ((t (/ j nt)))
              (when (>= j 1)
                (let ((newvalue (+ (- (psir 0)
                                      (* c1 (- (psir 1)
                                               (psir (- nx 1)))))
                                   (* c2 (+ (- (psir 1)
                                               (* 2 (psir 0)))
                                            (psir (- nx 1)))))))
                  (array-set! psinew newvalue 0))
                (let loop ((i 1))
                  (when (< i (- nx 1))
                    (let ((newvalue (+ (- (psir i)
                                          (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
                                       (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
                                                (psir (- i 1)))))))
                      (array-set! psinew newvalue i))
                    (loop (+ i 1))))
                (let ((newvalue (+ (- (psir (- nx 1))
                                      (* c1 (- (psir 0) (psir (- nx 2)))))
                                   (* c2 (+ (- (psir 0) (* 2 (psir (- nx 1))))
                                            (psir (- nx 2)))))))
                  (array-set! psinew newvalue (- nx 1))))
              (let loop ((i 0))
                (when (< i nx)
                  (array-set! psi (array-ref psinew i) i)
                  (array-set! X (list-ref x i) i j)
                  (array-set! Y t i j)
                  (array-set! Z (psir i) i j)
                  (loop (+ i 1)))))
            (loop (+ j 1)))))
      (list X Y Z)))))

(define (main args)
  (display "Hello World!\n")
  (let ((res (fd1d-advection-lax-wendroff)))
    (let ((X (list-ref res 0))
          (Y (list-ref res 1))
          (Z (list-ref res 2)))
      ; (display Z)
      (let ((port (open-output-pipe "python")))
        (format port "import pylab as pl\n")
        (format port "y = [[float(j) for j in i] for i in eval('~A'[2:].replace(' ',', '))]\n" Z)
        (format port "pl.imshow(y)\n")
        (format port "pl.colorbar()\n")
        (format port "pl.show()\n")
        (close-pipe port))))
  (newline))


-- 
Unpolitisch sein
heißt politisch sein
ohne es zu merken

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 298 bytes --]

Re: I wrote fluid advection code: How to make this more elegant?

[Panicz Maciej Godek]

[-- Attachment #1: Type: text/plain, Size: 6568 bytes --]

Hi,
although I cannot be of immediate help with the topic, because I don't know
anything about advection, I think that the main problem with your code is
that it is imperative. While Python's stylistics handles the imperative
code quite nicely, it looks rather strange in Scheme.

From my experience, the proper Scheme solution should resemble the
mathematical formulation of a problem (except that it should be more
descriptive), rather than a list of steps for solving it.

Also, it should be more readable to use patern matching instead of list
indexing, so most likely your expression

(let ((newvalue (+ (- (psir i)
                          (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
                       (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
                                (psir (- i 1)))))))
      (array-set! psinew newvalue i))

should look like this

(map (lambda (prev this next)
          (- this
             (* c1 (- next prev))
             (* (- c2) (+ next (* -2 this) prev))))
       `(0 ,@(drop psir 1))
       psir
       `(,@(drop-right psir 1) 0))

While this may also look slightly difficult to read (and write), this isn't
solely because of the expression's structure, but because the factors of
the expression have no name, and therefore the source code doesn't explain
their role (this is the problem of the Python code either, but Python
doesn't prompt to fix that)

HTH

PS I think that this subject fits better to guile-user


2016-01-23 11:00 GMT+01:00 Arne Babenhauserheide <arne_bab@web.de>:

> Hi,
>
> I just recreated a fluid advection exercise in Guile Scheme and I’m not
> quite happy with its readability. Can you help me improve it?
>
> My main gripe is that the math does not look instantly accessible.
>
> The original version was in Python:
>
>     psi[i] - c1*(psi[i+1] - psi[i-1]) + c2*(psi[i+1] - 2.0*psi[i] +
> psi[i-1])
>
> My port to Scheme looks like this:
>
>     (let ((newvalue (+ (- (psir i)
>                           (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
>                        (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
>                                 (psir (- i 1)))))))
>       (array-set! psinew newvalue i))
>
>
> Liebe Grüße,
> Arne
>
> Here’s the full code:
>
> #!/bin/sh
> # -*- scheme -*-
> exec guile -e '(@@ (advection) main)' -s "$0" "$@"
> !#
>
> ; Copyright (c) 2015 John Burkardt (original Python), 2016 Corinna
> ; Hoose (adaption) and 2016 Arne Babenhauserheide (pep8 + Scheme
> ; version).
>
> ; License: LGPL, built on the Python version from 2015 John Burkardt
> ; and Corinna Hoose. License LGPL.
>
> (define-module (advection)
>   #:use-module (ice-9 optargs) ; define*
>   #:use-module (srfi srfi-1) ; iota
>   #:use-module (ice-9 format)
>   #:use-module (ice-9 popen))
>
>
> (define* (fd1d-advection-lax-wendroff #:key (nx 101) (nt 1000) (c 1))
>   (let* ((dx (/ 1 (- nx 1)))
>          (x (iota nx 0 (/ 1 nx)))
>          (dt (/ 1 nt))
>          (c1 (* #e0.5 (* c (/ dt dx))))
>          (c2 (* 0.5 (expt (* c (/ dt dx)) 2))))
>     (format #t "CFL condition: dt (~g) ≤ (~g) dx/c\n" dt (/ dx c))
>     (let ((psi (make-array 0 nx))
>           (X (make-array 0 nx (+ nt 1)))
>           (Y (make-array 0 nx (+ nt 1)))
>           (Z (make-array 0 nx (+ nt 1))))
>       (let ((psinew (let ((pn (make-array 0 nx)))
>                       (let loop ((i 0))
>                         (cond ((= i nx)
>                                pn)
>                               (else
>                                (let ((xi (list-ref x i)))
>                                  (when (and (<= 0.4 xi) (<= xi 0.6))
>                                    (array-set! pn
>                                                (* (expt (- (* 10 xi) 4) 2)
>                                                   (expt (- 6 (* 10 xi)) 2))
>                                                i))
>                                  (loop (+ 1 i)))))))))
>         (define (psir i) (array-ref psi i))
>         (let loop ((j 0))
>           (cond
>            ((> j nt) #t) ; done
>            (else
>             (let ((t (/ j nt)))
>               (when (>= j 1)
>                 (let ((newvalue (+ (- (psir 0)
>                                       (* c1 (- (psir 1)
>                                                (psir (- nx 1)))))
>                                    (* c2 (+ (- (psir 1)
>                                                (* 2 (psir 0)))
>                                             (psir (- nx 1)))))))
>                   (array-set! psinew newvalue 0))
>                 (let loop ((i 1))
>                   (when (< i (- nx 1))
>                     (let ((newvalue (+ (- (psir i)
>                                           (* c1 (- (psir (+ i 1)) (psir (-
> i 1)))))
>                                        (* c2 (+ (- (psir (+ i 1)) (* 2
> (psir i)))
>                                                 (psir (- i 1)))))))
>                       (array-set! psinew newvalue i))
>                     (loop (+ i 1))))
>                 (let ((newvalue (+ (- (psir (- nx 1))
>                                       (* c1 (- (psir 0) (psir (- nx 2)))))
>                                    (* c2 (+ (- (psir 0) (* 2 (psir (- nx
> 1))))
>                                             (psir (- nx 2)))))))
>                   (array-set! psinew newvalue (- nx 1))))
>               (let loop ((i 0))
>                 (when (< i nx)
>                   (array-set! psi (array-ref psinew i) i)
>                   (array-set! X (list-ref x i) i j)
>                   (array-set! Y t i j)
>                   (array-set! Z (psir i) i j)
>                   (loop (+ i 1)))))
>             (loop (+ j 1)))))
>       (list X Y Z)))))
>
> (define (main args)
>   (display "Hello World!\n")
>   (let ((res (fd1d-advection-lax-wendroff)))
>     (let ((X (list-ref res 0))
>           (Y (list-ref res 1))
>           (Z (list-ref res 2)))
>       ; (display Z)
>       (let ((port (open-output-pipe "python")))
>         (format port "import pylab as pl\n")
>         (format port "y = [[float(j) for j in i] for i in
> eval('~A'[2:].replace(' ',', '))]\n" Z)
>         (format port "pl.imshow(y)\n")
>         (format port "pl.colorbar()\n")
>         (format port "pl.show()\n")
>         (close-pipe port))))
>   (newline))
>
>
> --
> Unpolitisch sein
> heißt politisch sein
> ohne es zu merken
>

[-- Attachment #2: Type: text/html, Size: 8430 bytes --]

Re: I wrote fluid advection code: How to make this more elegant?

[Arne Babenhauserheide]

[-- Attachment #1: Type: text/plain, Size: 2363 bytes --]


Panicz Maciej Godek writes:
> although I cannot be of immediate help with the topic, because I don't know
> anything about advection, I think that the main problem with your code is
> that it is imperative. While Python's stylistics handles the imperative
> code quite nicely, it looks rather strange in Scheme.
>
> From my experience, the proper Scheme solution should resemble the
> mathematical formulation of a problem (except that it should be more
> descriptive), rather than a list of steps for solving it.

That’s a good point — I should have given the mathematical formulation
along with the code. The advection equation this implements is:

(φ_j,n+1 - φ_j,n)/Δt
  + c(φ_j+1,n - φ_j-1,n)/2Δx
  - (c² Δt/Δx) * (φ_j+1,n - 2φ_j,n + φ_j-1,n) / 2Δx
  = 0

[sadly this isn’t easy to read, either, but it’s the math which needs implementation]

> Also, it should be more readable to use pattern matching instead of list
> indexing, so most likely your expression
>
> (let ((newvalue (+ (- (psir i)
>                           (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
>                        (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
>                                 (psir (- i 1)))))))
>       (array-set! psinew newvalue i))
>
> should look like this
>
> (map (lambda (prev this next)
>           (- this
>              (* c1 (- next prev))
>              (* (- c2) (+ next (* -2 this) prev))))
>        `(0 ,@(drop psir 1))
>        psir
>        `(,@(drop-right psir 1) 0))
>
> While this may also look slightly difficult to read (and write),

I think it already helps a lot, that this is shorter (prev this next
clearly wins against (+/- i 1)).

> this isn't solely because of the expression's structure, but because
> the factors of the expression have no name, and therefore the source
> code doesn't explain their role (this is the problem of the Python
> code either, but Python doesn't prompt to fix that)

that part is already in the math, but using variables with better names
sounds like an improvement.

Thank you!

> PS I think that this subject fits better to guile-user

You’re right — sorry. I’m still spoiled by projects where exactly one
mailing list is active…

Best wishes,
Arne
-- 
Unpolitisch sein
heißt politisch sein
ohne es zu merken

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 298 bytes --]

Re: I wrote fluid advection code: How to make this more elegant?

[Taylan Ulrich Bayırlı/Kammer]

Arne Babenhauserheide <arne_bab@web.de> writes:

> Hi,
>
> I just recreated a fluid advection exercise in Guile Scheme and I’m not
> quite happy with its readability. Can you help me improve it?
>
> My main gripe is that the math does not look instantly accessible.
>
> The original version was in Python:
>
>     psi[i] - c1*(psi[i+1] - psi[i-1]) + c2*(psi[i+1] - 2.0*psi[i] + psi[i-1])
>
> My port to Scheme looks like this:
>
>     (let ((newvalue (+ (- (psir i)
>                           (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
>                        (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
>                                 (psir (- i 1)))))))
>       (array-set! psinew newvalue i))

Guile supports SRFI-105, so that could be:

    {{psi[i] - {c1 * {psi[{i + 1}] - psi[{i - 1}]}}} + {c2 * {{psi[{i + 1}] - {2 * psi[i]}} + psi[{i - 1}]}}}

which is less readable than the Python version because there's no
first-hand support for operator precedence so we {} group all binary
operations, plus we need to use {} within [], plus we must separate
operators with whitespace, but maybe it's acceptable.

You can put "#!curly-infix" at the start of your Scheme file to enable
SRFI-105 syntax.

Note that all those 'psi[x]' expressions will expand to

    ($bracket-apply$ psi x)

and you can implement a macro of that name to turn that into whatever is
suitable at compile-time.  (If performance is not an issue, SRFI-123
provides a run-time polymorphic procedure for $bracket-access$.)

With a smart definition of $bracket-apply$, you could also cut down
those psi[{i + 1}] to psi[i + 1].  The latter will expand to

    ($bracket-apply$ psi i + 1)

which could be accommodated for in a special-purpose definition of
$bracket-apply$ such as:

    (syntax-rules ()
      ((_ obj index)
       (obj index))
      ((_ obj x operator y)
       (obj (operator x y))))

That would *not* support e.g. psi[i + j - 1], which would instead need
to be written e.g. psi[{i + j} - 1], i.e. we only save one pair of {},
but maybe that's good enough.

By the way, e.g. {x - y + z} will turn into ($nfx$ x - y + z), and maybe
there's a definition for $nfx$ out in the wild (or you can create one)
that does operator precedence.  (That could then also be used in
$bracket-apply$.)  So optimally, the whole thing could become:

    {psi[i] - c1 * {psi[i + 1] - psi[i - 1]} + c2 * {psi[i + 1] - 2 * psi[i] + psi[i - 1]}}

Taylan

Re: I wrote fluid advection code: How to make this more elegant?

[Nala Ginrut]

The math formula may take advantage of sweet expression which is
contained in the current Guile reader, IIRC.
One may give it a try?

On Sat, 2016-01-23 at 13:42 +0100, Taylan Ulrich Bayırlı/Kammer wrote:
> Arne Babenhauserheide <arne_bab@web.de> writes:
> 
> > Hi,
> >
> > I just recreated a fluid advection exercise in Guile Scheme and I’m not
> > quite happy with its readability. Can you help me improve it?
> >
> > My main gripe is that the math does not look instantly accessible.
> >
> > The original version was in Python:
> >
> >     psi[i] - c1*(psi[i+1] - psi[i-1]) + c2*(psi[i+1] - 2.0*psi[i] + psi[i-1])
> >
> > My port to Scheme looks like this:
> >
> >     (let ((newvalue (+ (- (psir i)
> >                           (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
> >                        (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
> >                                 (psir (- i 1)))))))
> >       (array-set! psinew newvalue i))
> 
> Guile supports SRFI-105, so that could be:
> 
>     {{psi[i] - {c1 * {psi[{i + 1}] - psi[{i - 1}]}}} + {c2 * {{psi[{i + 1}] - {2 * psi[i]}} + psi[{i - 1}]}}}
> 
> which is less readable than the Python version because there's no
> first-hand support for operator precedence so we {} group all binary
> operations, plus we need to use {} within [], plus we must separate
> operators with whitespace, but maybe it's acceptable.
> 
> You can put "#!curly-infix" at the start of your Scheme file to enable
> SRFI-105 syntax.
> 
> Note that all those 'psi[x]' expressions will expand to
> 
>     ($bracket-apply$ psi x)
> 
> and you can implement a macro of that name to turn that into whatever is
> suitable at compile-time.  (If performance is not an issue, SRFI-123
> provides a run-time polymorphic procedure for $bracket-access$.)
> 
> With a smart definition of $bracket-apply$, you could also cut down
> those psi[{i + 1}] to psi[i + 1].  The latter will expand to
> 
>     ($bracket-apply$ psi i + 1)
> 
> which could be accommodated for in a special-purpose definition of
> $bracket-apply$ such as:
> 
>     (syntax-rules ()
>       ((_ obj index)
>        (obj index))
>       ((_ obj x operator y)
>        (obj (operator x y))))
> 
> That would *not* support e.g. psi[i + j - 1], which would instead need
> to be written e.g. psi[{i + j} - 1], i.e. we only save one pair of {},
> but maybe that's good enough.
> 
> By the way, e.g. {x - y + z} will turn into ($nfx$ x - y + z), and maybe
> there's a definition for $nfx$ out in the wild (or you can create one)
> that does operator precedence.  (That could then also be used in
> $bracket-apply$.)  So optimally, the whole thing could become:
> 
>     {psi[i] - c1 * {psi[i + 1] - psi[i - 1]} + c2 * {psi[i + 1] - 2 * psi[i] + psi[i - 1]}}
> 
> Taylan
> 

Re: I wrote fluid advection code: How to make this more elegant?

[Arne Babenhauserheide]

[-- Attachment #1: Type: text/plain, Size: 2347 bytes --]


Taylan Ulrich Bayırlı/Kammer writes:
> Arne Babenhauserheide <arne_bab@web.de> writes:
>> The original version was in Python:
>>
>>     psi[i] - c1*(psi[i+1] - psi[i-1]) + c2*(psi[i+1] - 2.0*psi[i] + psi[i-1])
>>
>> My port to Scheme looks like this:
>>
>>     (let ((newvalue (+ (- (psir i)
>>                           (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
>>                        (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
>>                                 (psir (- i 1)))))))
>>       (array-set! psinew newvalue i))
>
> Guile supports SRFI-105, so that could be:
>
>     {{psi[i] - {c1 * {psi[{i + 1}] - psi[{i - 1}]}}} + {c2 * {{psi[{i + 1}] - {2 * psi[i]}} + psi[{i - 1}]}}}

That’s already pretty close — I wonder why I didn’t think of the psi[i]
form.

I think a + around the equation would actually help here:

    (+ psi[i]
       (* -1 c1 {psi[{i + 1}] - psi[{i - 1}]})
       (* c2 {{psi[{i + 1}] - {2 * psi[i]}} + psi[{i - 1}]}))

Though neoteric expressions combined with curly infix make this even
easier: p{i + 1} → (p (+ i 1))

(though this did not work for me in the REPL right now — did I miss something?)

So the function psir could be used to have elegant access to elements:

    (+ psi[i]
       (* -1 c1 {psi[{i + 1}] - psi[{i - 1}]})
       (* c2 {{psi[{i + 1}] - {2 * psi[i]}} + psi[{i - 1}]}))

>     {psi[i] - c1 * {psi[i + 1] - psi[i - 1]} + c2 * {psi[i + 1] - 2 * psi[i] + psi[i - 1]}}

That looks roughly as readable as the Python version. With the + around
I think it becomes better:

    (+ psi[i]
       (* -1 c1 {psi[i + 1] - psi[i - 1]})
       (* c2 {{psi[i + 1] - {2 * psi[i]}} + psi[i - 1]}))

And I don’t think it’s worse than the Python version, because in Python
I also need to first import numpy which changes some semantics (the []
operator changes and arrays have some special properties which can
sometimes bite).

Besides: when working on the Python, I was tempted to bracket the
expression to make it easier to split the lines and make the boundary
between sub-expressions clearer:

    (psi[i]
      - c1*(psi[i+1] - psi[i-1])
      + c2*(psi[i+1] - 2.0*psi[i] + psi[i-1]))

Best wishes,
Arne

PS: CCed to David A. Wheeler who wrote SRFI-105.
-- 
Unpolitisch sein
heißt politisch sein
ohne es zu merken

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 298 bytes --]

Re: I wrote fluid advection code: How to make this more elegant?

[Taylan Ulrich Bayırlı/Kammer]

Arne Babenhauserheide <arne_bab@web.de> writes:

> Taylan Ulrich Bayırlı/Kammer writes:
>> Arne Babenhauserheide <arne_bab@web.de> writes:
>>> The original version was in Python:
>>>
>>>     psi[i] - c1*(psi[i+1] - psi[i-1]) + c2*(psi[i+1] - 2.0*psi[i] + psi[i-1])
>>>
>>> My port to Scheme looks like this:
>>>
>>>     (let ((newvalue (+ (- (psir i)
>>>                           (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
>>>                        (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
>>>                                 (psir (- i 1)))))))
>>>       (array-set! psinew newvalue i))
>>
>> Guile supports SRFI-105, so that could be:
>>
>>     {{psi[i] - {c1 * {psi[{i + 1}] - psi[{i - 1}]}}} + {c2 * {{psi[{i + 1}] - {2 * psi[i]}} + psi[{i - 1}]}}}
>
> That’s already pretty close — I wonder why I didn’t think of the psi[i]
> form.
>
> I think a + around the equation would actually help here:
>
>     (+ psi[i]
>        (* -1 c1 {psi[{i + 1}] - psi[{i - 1}]})
>        (* c2 {{psi[{i + 1}] - {2 * psi[i]}} + psi[{i - 1}]}))
>
> Though neoteric expressions combined with curly infix make this even
> easier: p{i + 1} → (p (+ i 1))

Good call, I had forgotten that's included in SRFI 105.  (Thought it was
just for x[y].)

> (though this did not work for me in the REPL right now — did I miss something?)

Any SRFI 105 syntax must appear within {}, so that would have to be
e.g. {p{i + 1}}, although if it appears within a larger {} block then it
won't have that annoyance.

> So the function psir could be used to have elegant access to elements:
>
>     (+ psi[i]
>        (* -1 c1 {psi[{i + 1}] - psi[{i - 1}]})
>        (* c2 {{psi[{i + 1}] - {2 * psi[i]}} + psi[{i - 1}]}))
>
>>     {psi[i] - c1 * {psi[i + 1] - psi[i - 1]} + c2 * {psi[i + 1] - 2 * psi[i] + psi[i - 1]}}
>
> That looks roughly as readable as the Python version. With the + around
> I think it becomes better:
>
>     (+ psi[i]
>        (* -1 c1 {psi[i + 1] - psi[i - 1]})
>        (* c2 {{psi[i + 1] - {2 * psi[i]}} + psi[i - 1]}))

I agree mixing in some prefix notation actually makes things more
readable here.

Taylan