Evaluating Polynomials at Grid Points

I am continuing my quest on B-splines. The code below builds a 5x5 matrix out of B-splines, using the BSplineBasis[] routine.

I now want to evaluate the polynomials that are stored in each matrix element. The problem that I am facing is, however, that each matrix element corresponds to a different position on my xy-grid.

For instance, B[[2,2]] corresponds to the grid-pint x=-1, y = -1 . B[[2,3]] corresponds to x=-1, y = 0 and so forth.

At the moment I do x = smth and then evaluate the polynomial manually. This works for a small system, but it is not optimal for larger problems.

So my question is, is there some notation that does something like B[[2,3]](put out the value of this matrix element at (x = smth,y = smth))? I tried it with Evaluate[], but didn't get it to work.

knots = {-1, -1, -1 , -1, 0, 1, 1, 1, 1};
B1 = Table[
D[BSplineBasis[{3, knots}, i, x], {x, 2}] BSplineBasis[{3, knots},
j, y], {i, 0, 4}, {j, 0, 4}];
B2 = Table[
BSplineBasis[{3, knots}, i, x] D[
BSplineBasis[{3, knots}, j, y], {y, 2}], {i, 0, 4}, {j, 0, 4}];
B = B1 + B2;
(*set boundary conditions*)
B[[1]] = {0, 0, 0, 0, 0};
B[[5]] = {0, 0, 0, 0, 0};
B[[All, 5]] = {0, 0, 0, 0, 0};
B[[2, 1]] = BSplineBasis[{3, knots}, 2, x];
B[[3, 1]] = BSplineBasis[{3, knots}, 3, x];
B[[4, 1]] = BSplineBasis[{3, knots}, 4, x];


Here is a plot of what I'm doing, just for nurturing the eye.

Plot3D[Evaluate[B], {x, -1, 1}, {y, -1, 1}, Mesh -> None,
PlotStyle -> Directive[Opacity[0.8]], PlotRange -> {-2, 2},
ColorFunction -> Function[{x, y, z}, Hue[z]]]


-
How do we know that B[[2,3]] corresponds to the point (-1,0)? It seems that you want to evaluate B[[2,3]] at that point. Is that correct? Can you just run B[[2,3]]/.{x->-1,y->0}? –  Jonathan Shock May 1 '13 at 10:16
@JonathanShock You can not know it without giving some meaning to it. In this case the meaning is given by the physical problem that I am trying to solve. In itself, this choice is completely arbitrary. Thanks! What you posted is exactly what I needed. If you want to put it in as an answer, I can upvote. –  lomppi May 1 '13 at 10:18
thanks. I've done so. –  Jonathan Shock May 1 '13 at 10:23
no, thank you!!! –  lomppi May 1 '13 at 10:28
I get a pile of Power::infy: Infinite expression 1/0 encountered. >> errors when I try to run your code. Can you correct that? –  Mr.Wizard May 1 '13 at 10:33

You can perform a replace:

B[[2,3]]/.{x->-1,y->0}


If you have a table of x and y values (call it xyvals) for given matrix elements then you could do:

Table[B[[m,n]]/.{x->xyvals[[m,n,1]],y->xyvals[[m,n,2]]},{m,Length[B]},{n,Length[B[[1]]]}]

-
I tried the Table thing you posted. So I did: xyvals = Table[{0, 0}, {i, 1, 5}, {j, 1, 5}] just trying to set it up. However, I get error messages for simple cases such as: xyvals[[2, 2, 1]]. I know this is simple stuff, please bear with me, I'm coming to Mathematica from Python. –  lomppi May 1 '13 at 11:58
In your example you don't have a B[[2,2]]. Could this be the problem? If you evaluate xyvals[[2,2,1]] does it give you an error? Can you state exactly what you're evaluating which gives the problem? –  Jonathan Shock May 1 '13 at 12:04
no, it's much simpler than that. I just don't know how to use Table[] ;-). but I'm learning. trying out examples from the Ref to get a better understanding right now –  lomppi May 1 '13 at 12:13
I think unless you have a specific example Ref is going to be your best friend. –  Jonathan Shock May 1 '13 at 12:24

One approach is to turn your expression into functions. I give a few extra variations just to show what can be done and which you might find useful to learn.

Here's the whole thing as a function (to be used, for example, in your Plot3D):

Bfn = Evaluate[B /. {x -> #1, y -> #2}] &;

Bfn[x, y] == B
(* True *)


Component functions can be extracted from Bfn,

BPartFn[i_, j_] := Evaluate@Bfn[[1, i, j]] &


or from B,

BPartFn[i_, j_] := Evaluate[B[[i, j]] /. {x -> #1, y -> #2}] &


Comparisons:

BPartFn[2, 3][x, y] == B[[2, 3]]
(* True *)

BPartFn[2, 3][-1, 0]
(* -9/2 *)

B[[2, 3]] /. {x -> -1, y -> 0}
(* -9/2 *)


If you need to evaluate them at many points, you might wish to make each function Listable:

BPartFn[i_, j_] := Function[Null, Evaluate@Bfn[[1, i, j]], Listable]


Then the function will automatically thread over lists of coordinates with a call of the form

BPartFn[2, 3][ {x1, x2, ...}, {y1, y2, ...} ]


Here are 10 points:

xy = RandomReal[{0, 1}, {10, 2}]
(* {{0.0214747, 0.362779}, {0.00453566, 0.0420441}, {0.126681, 0.467454},
{0.422858, 0.155077}, {0.482303, 0.920426}, {0.422345, 0.70505}, {0.572743, 0.121306},
{0.912936, 0.601003}, {0.39634, 0.437055}, {0.542907, 0.867795}} *)


Comparison of outputs:

BPartFn[2, 3] @@ Transpose[xy]
(* {0.321359, 0.0651208, 0.326909, 0.305393, 0.0944842, 0.150124,
0.263136, 0.0229906, 0.248124, 0.0690792} *)

B[[2, 3]] /. (Thread[{x, y} -> #] & /@ xy)
(* {0.321359, 0.0651208, 0.326909, 0.305393, 0.0944842, 0.150124,
0.263136, 0.0229906, 0.248124, 0.0690792} *)


If you're going to be making a lot of calls to BPartFn, you can get a little improvement in performance by memoizing:

    BPartFn[i_, j_] := BPartFn[i_, j_] = Evaluate@Bfn[[1, i, j]] &


(It can be done with the other definitions as well.)

-
thanks for this. this is quite impressive. do you know of any tutorial that teaches this? I tried wolfram, but couldn't find anything comprehensive. –  lomppi May 1 '13 at 19:23
@sebastian The Core Language tutorial might be a good place to look (there are a couple of chapters on functions). This question aims to be an index to all sorts of resources. I think many of us learn to select the function or operator and use the Help > Find Selected Function menu command, when we see something we don't understand; and then root around. –  Michael E2 May 1 '13 at 20:36