Plotting a lambda that returns a (math/matrix) vector

Questions & Answers
,
1

I'm writing a text where I'm needing to plot polynomial bases and operations on them. The bases are often represented in a vector form, e.g.,

\mathbf{N}^p(\xi) = [ N_0^p(\xi), N_1^p(\xi), \dots, N_p^p(\xi) ]^T

Where N_i^p(\xi) may be a power, Bernstein, Lagrange, Legendre, Chebyshev, etc. basis function

I then want to be able to do things like plot \mathbf{N}^p(\xi) for \xi \in [\xi_{\min}, \xi_{\max}] , or multiply \mathbf{N}^p(\xi) by a matrix and plot the result (\mathbf{S}^p(\xi) below):

\mathbf{S}^p(\xi) = \left[\mathbf{C}\right]^{n \times (p+1) } \left( \mathbf{N}^p(\xi) \right)

Code for assembling a math/matrix of a Bernstein basis:

#lang racket/base
(require math)
(require plot)

(define (evalBernsteinBasis1D degree basis-idx variate)
  (define term1 (binomial degree basis-idx) )
  (define term2 (expt variate basis-idx) )
  (define term3 (expt (- 1 variate) (- degree basis-idx) ) )
  (* term1 term2 term3)
)

(define (eval-bernstein-basis-1D-vector degree variate ) 
  (define basis-list (for/list ( [i (in-range (+ degree 1 ))] ) (evalBernsteinBasis1D degree (+ i 0) variate) ) )
  (list->matrix (+ degree 1) 1 basis-list )
)

Working, but verbose plotting of vector

(parameterize ( [plot-width    800]
                [plot-height   400]
                [plot-x-label  #f]
                [plot-y-label  #f]
                [plot-pen-color-map `tab10])
  (plot
   (list
    (function (lambda (x) (matrix-ref (eval-bernstein-basis-1D-vector 2 x ) 0 0) ) 0 1 #:y-min 0 #:y-max 1 #:color 0 #:width 4)
    (function (lambda (x) (matrix-ref (eval-bernstein-basis-1D-vector 2 x ) 1 0) ) 0 1 #:y-min 0 #:y-max 1 #:color 0 #:width 4)
    (function (lambda (x) (matrix-ref (eval-bernstein-basis-1D-vector 2 x ) 2 0) ) 0 1 #:y-min 0 #:y-max 1 #:color 0 #:width 4)
    )
   )
)

Results in:

plot
plot800×400 35.5 KB

But I'd like to know if I can be more efficient with this plotting - a single function plot and no (matrix-ref _ 0 0)

For example,

An attempt at a more efficient plot

(parameterize ( [plot-width    800]
                [plot-height   400]
                [plot-x-label  #f]
                [plot-y-label  #f]
                [plot-pen-color-map `tab10])
  (plot
   (function (lambda (x) (eval-bernstein-basis-1D-vector 2 x ) ) 0 1 #:y-min 0 #:y-max 1 #:color 0 #:width 4)
  )
)

But this results in an empty plot.

2

The function renderer needs a function which takes a number and returns a single number, so it won't work to return a 3-element matrix and expect it to plot 3 lines -- I am not sure if any plot library can do it like that. On the other hand, it seems that eval-bernstein-basis-1D-vector just calls evalBernsteinBasis1D 3 times to create a matrix, than you just extract one element from the matrix, so it might be simpler to just create the 3 plot renderers directly:

(define (make-bernstein-renderer degree x-min x-max)
  (for/list ([basis-idx (in-range (add1 degree))])
    (function
     (lambda (x)
       (evalBernsteinBasis1D degree basis-idx x))
     x-min x-max #:color 0)))

(plot (make-bernstein-renderer 2 0 1))

Another option is to do what numerical plotting libraries such as Matplotlib or Octave do and generate the entire data set to be plotted:

;; Evaluate a polynomial, POL to DEGREE between VARIATE-MIN and VARIATE-MAX
;; generating (N + 1) samples.
(define (eval-pol pol degree variate-min variate-max n)
  (define step (/ (- variate-max variate-min) n))
  (build-matrix
   (add1 degree)                        ; rows
   (add1 n)                             ; columns
   (lambda (row col)
     (define variate (+ variate-min (* col step)))
     (pol degree row variate))))

For example, (eval-pol evalBernsteinBasis1D 2 0 1 10) will create a matrix of 11 samples between 0 and 1:

(array
 #[#[1 81/100 16/25 49/100 9/25 1/4 4/25 9/100 1/25 1/100 0]
   #[0 9/50 8/25 21/50 12/25 1/2 12/25 21/50 8/25 9/50 0]
   #[0 1/100 1/25 9/100 4/25 1/4 9/25 49/100 16/25 81/100 1]])

You can get obtain the values from an individual basis using (array-axis-ref a 0 basis-index), which is a "horizontal slice", for example:

> (array-axis-ref a 0 1)
(array #[0 9/50 8/25 21/50 12/25 1/2 12/25 21/50 8/25 9/50 0])

... or you can extract the entry for a variate (a vertical slice in the array), for example:

> (array-axis-ref a 1 1) ; get the row for the second variate
(array #[81/100 9/50 1/100])

With some helpers you can plot this data numericaly, just as you would do with numpy and Matplotlib:

;; Plot a lines using X and Y arrays, similar to how Matplot lib does line plots
(define (line-plot-xy x y #:color (c 0))
  (lines
   (for/list ([xval (in-array x)]
              [yval (in-array y)])
     (vector xval yval))
   #:color c))

;; Generate the X values between 0 and 1 (11 samples):
(define x (linspace 0 1 10))
(define degrees 2)
;; Evaluate Bernstein with degree 2 between 0 and 1 with 11 samples:
(define b (eval-pol evalBernsteinBasis1D degrees 0 1 10))

;; plot each individual row in the resulting data set:
(plot
 (for/list ([basis-index (in-range (add1 degrees))])
   (line-plot-xy x (array-axis-ref b 0 basis-index))))

Hope this helps,
Alex.

3

Alex has adressed how to use plot.

In case you want to plot several polynomials given as coordinates in a shared base, it makes sense to compute the common base first, and then use matrix-dot between the coordinates and the base.

4

My two "inspirations" were Matlab's symbolic fplot routine:

clear
x = sym("x");
B = [x; x^2; x^3];
fplot(B, [0 1])

and the more traditional plot routine with an anonymous function, evaluated to an array:

clear
B = @(x) [ x; x.^2; x.^3 ]
x = linspace(0, 1, 1e3);
plot(x,B(x))

Thanks for your two examples of different approaches - I'll explore both examples in more detail to decide which to go with.

1 Like