Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]]]

enter image description here

share|improve this question
1  
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
show 2 more comments

2 Answers

up vote 3 down vote accepted

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]]]}]
share|improve this answer
    
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
add comment

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.)

share|improve this answer
    
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
1  
@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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.