Enzo,

I think you found a bug in Mathematica's BSplineBasis. It should be putting a conditional on the terms of the derivative. I would report this at Wolfram Tech Support or better yet, you should call them and see if they can suggest a good workaround until it is fixed. I added a few lines to your file showing the problematic expression. The derivative has two terms and the first term will go to zero in the denominator when n->0. The second term has the answer that agrees with your answer when you substitute in zero for n and then take the derivative. Note that the derivative expression also goes to infinity when n > 5. In this case the second term of the derivative has zero in its denominator.

As for a workaround:

D[BSplineBasis[{d, U}, n, x], {x,1}] appears to work properly if the n value is numeric before you take the derivative. I suggest you write your code so that you make sure it happens that way. One possibility is to use operator overloading to compute all the derivative expressions as functions of x and then evaluate them.

If you define:

bSplineDeriv[n_] := 
 bSplineDeriv[n] = D[BSplineBasis[{p, U}, n, x], {x, 1}]

As you iterate through your values of n, the symbolic expressions for the derivative will be computed. The nice part of this definition is that it will only compute the derivative for each n value once. See the enclosed nb file.

I hope this helps. Maybe someone else has another workaround suggestion...

Attachments:

Enzo,

Next, try to do the procedure exactly as you would do in Matlab but using Mathematica to handle the equations. Enter the PDE equations, fit splines, take derivatives and construct large matrices -- Construct functions that you use for each step but don't defer the computation. For example, Don't create a formula to construct spines and use it later. Define the formula for splines and create a matrix of splines (I may be showing my lack of understanding of your algorithm here -- but you get the idea). After you have an example that works, then you can create functions that do the steps in that order without delaying the calculations until the end (which causes repetition). If you have done this in Matlab it should be relatively straightforward in Mathematica. If you post some Matlab code I might be able to suggest more details but essentially you should let Mathematica handle the large expressions and expand everything out into symbolic or numerical values early -- at the same point in the procedure that you would do so in Matlab. It is better to consume memory than to repeat your calculations so it is better to iterate and create a large matrix of expressions and then take derivatives of the matrix and later turn it into numbers than to defer the calculation and do the derivatives and expression handling while you are creating the large numerical matrix. (the later may seem more "elegant" but it will be slower).

I hope that helps.

Regards

Thank you Neil for your comments! Yes, I see the definition `bspU[i_]:=bspU[i] helps! And I agree with you: I think that a lot of time can be saved by avoiding Mathematica to compute the derivative everytime with command D. The derivative of a bspline can be defined recursively form the bspline itself. Unfortunately, I haven't found it between the Mathematica built-in functions so I will try to define it by myself hopig that it will be faster than computing it via D.

As regards Sander's version of my code, it is faster. The only thing is that I need to let Mathematica to collect terms since in the real applications the original differential equations is rather lengthy and qiute complex.

Many thanks!

Enzo

Hello Neil,

Thank you very much for your helpful and instructive reply to my post. The example that you kindly provided helps a lot to understand the issue.
Actually, I also had the feeling that something was wrong since all the functions remain unevaluated until the very end and they are called (unnecessarily) many times, but I did not know how to face it. I will try to reformulate my code taking inspiration from the approach you suggest in your example. I will let you know.

Many thanks!

Enzo

Again, tell us what K is. What does it represent? I know K is a matrix, but what do the elements represent? Without knowing any background it is SUPER HARD to come up with a faster solution. Your code is also very hard to read because of all the functions that nearly don't do anything, very convoluted. I cleaned up your code a lot, which took me a lot of time to figure out what you are doing!

p=5;
q=5;
n=8;
m=8;
r=n+p+1;
s=m+q+1;
rint=r+1-2 (p+1);
sint=s+1-2 (p+1);
Uint=N@Subdivide[rint+1];
Vint=N@Subdivide[sint+1];
U=Join[ConstantArray[0.0,p],Uint,ConstantArray[1.0,p]];
V=Join[ConstantArray[0.0,q],Vint,ConstantArray[1.0,q]];
Uc=MovingAverage[U[[2;;-2]],p];
Vc=MovingAverage[V[[2;;-2]],p];

ClearAll[DoIt,D1w11V2,D2w112V2,bspU,bspV]

bspU[i_]:=bspU[i]=BSplineBasis[{p,U},i-1,u]
bspV[j_]:=bspV[j]=BSplineBasis[{q,V},j-1,v]

D1w11V2[{i_,j_}]:=D[bspU[i],u]bspV[j]
D2w112V2[{i_,j_}]:=D[bspU[i],u] D[bspV[j],v]

DoIt[rr_List,ij:{i_,j_}]:=Module[{eq},
    eq=D1w11V2[ij]+D2w112V2[ij];
    Table[eq/.{u->Uc[[r[[1]]]],v->Vc[[r[[2]]]]},{r,rr}]
]

ivec=Flatten[ConstantArray[Range[n+1],m+1]];
jvec=Flatten[Transpose[ConstantArray[Range[m+1],n+1]]];
ijvec={ivec,jvec}\[Transpose];

K=Transpose[Table[DoIt[ijvec,s],{s,ijvec}]];

Total[Abs[Chop[(K - Korig), 10^-10]], \[Infinity]] (* check if K and original K are the same within numerical error *)

This runs already a bit faster, but without knowing what you're doing, mathematically, or even in spoken words, I'm afraid it will be very tough to speed it up. What does a row represent in 'k'? what does a column represent? Also posting your Matlab code would be helpful maybe... Attach it to your post.

Hello Sander,

Thank you so much for your reply and for the time you spent in cleaning up my code.

Here the mathematical meaning of what I need to do:

I want to solve a differential equation (d.e.) through a numerical method which basically transforms the d.e. into an algebraic system in the standard form Kx = b. Once I find K, then I will just need to invert it and get the unknown vector x = K^-1 b. For the sake of simplicity, I have posted a very dummy example of differential equation. (Please, note that it is a meaningless example since neither the term b nor the boundary conditions are present). The d.e. of the dummy example reads as follows

w1,1 + w1,12 = 0 (1)

where w1(u,v) is a function defined from a square domain [0,1] x [0,1] to R. In the code that I posted D1w[1, 1] is w1,1, namely the partial derivative of w1 with respect to u, while D2w[1, 1, 2] is w1,12, namely the second derivative of w1 with respect to u and v. There are several ways to approximate w1. I have chosen to approximate it by means of BSpline basis functions. So I have

w1(u,v) = \Sum{i=0,n} \Sum{j=0,m} Ni(u)*Nj(v) W1(i,j) (2)

where W1(i,j) represent a (finite) set of (n+1)(m+1) control variables, and Ni and Nj are the i-th and j-th BSpline basis functions. Now, it is easy to take the derivatives of the approximate form of w1 as given in Eq. (2) and replace them into the original d.e. Eq (1). If one makes the substitution, which is what I try to do with `EqDisc[ic, jc, i, j] := Eq //. sub[ic, jc, i, j]; (see my code) one gets a discretized version of Eq (1), which represents one (approximated) equation in (n+1)(m+1) unknowns.

If now one evaluates this equation in (n+1)(m+1) points belonging to the domain [0,1]x[0,1] (let’s say a grid of points) one gets (n+1)(m+1) equations. A linear system is obtained and, after collecting terms (see command Coefficient), the matrix K of the system is found. So, each row of K represents the discretized equation evaluated in a specific point of the domain, while each element on a row represents the term appearing in the summations over i and j. Please note that a conversion from two-index to one-index is needed since the unknowns W1(i,j) are in the final system arranged in a vector and not in a matrix.

Sorry for the very informal explanation, but I hope this is enough to give you a general idea of what my final goal is.

This is a known technique which works well when I implement it in Matlab environment where I use a more standard programming approach. Now I would like to develop the same technique in Mathematica.

Many thanks, Enzo

Hi Enzo,

the codes takes 225 milliseconds on my machine, 1000x faster? 225 microseconds? It would be very nice to explain WHAT you are doing/what you want to accomplish. A literal code-translation is generally not a good way of doing things efficiently in Mathematica.

I was also wondering why you have:

 Flatten[ArrayReshape[ ... ]]

The ArrayReshape doesn't do anything, right? Also using Table to make copies is slower than when you use ConstantArray. In some occasions you use //. while a single replacement /. is probably enough right?

Uint = Range[0, 1, 1/(rint + 1)]; 
Vint = Range[0, 1, 1/(sint + 1)];

I think SubDivide is the function you want to use here, faster and more legible. Lastly, I see that the end result is part integers, part real numbers. Try to 'plug in' real values as soon as you can, this can speed up code by a very large amount, also having only reals (not mixed types) speeds things up a lot.

Again, tell us what you (want to) calculate, now it is difficult to guess because of all the functions that call each other...