Enzo,

I did not extensively study your approach but I think I see your problem.

Easy question first. To make a number a float just add the decimal point or a prime (above the tab key) -- instead of 1,2,3 you use 1.,2.,3. or 1.0,2.0,3.0 (because it looks better!), or 1`, 2`, 3` (The prime ties into the Precision)

When you switch from Matlab to Mathematica you need to change the way you think about solving a problem. In general, using Matlab, matrices and coefficients are computed as you proceed because it is a matrix-based language. Mathematica is primarily symbolic. The way you handled these calculations is you created functions. By creating Mathematica functions you delay the calculation of everything until you call it. Therefore many of your calculations are duplicated MANY times. You need to rework your approach a bit to make this faster. For a simple example (notebook attached):

Create a large equation

In[1]:= exp = 1000;

In[2]:= eq = (x + 3)^exp

Out[2]= (3 + x)^1000

Compute the derivative with the large equation (just to have something to compute)

In[3]:= eq2 = D[eq, x]/exp

Out[3]= (3 + x)^999

Create a function that takes the derivative -- and slow it down some more by Expanding first

In[4]:= computefun[f_, var_] := D[Expand[f], var]

Create two functions -- one that calls the compute function and one that takes the computed result and make a function

In[5]:= fun1[i_, var_] := Coefficient[computefun[eq, var], var, i]

In[6]:= fun2[i_] := Coefficient[eq2, x, i]

In[7]:= 

Slow -- it recomputes the result

In[8]:= Timing[mat = Table[fun1[i, x], {i, exp}];]

Out[8]= {2.20229, Null}

Faster -- it calculates it once using table

In[9]:= Timing[mat = Table[fun2[i], {i, exp}];]

Out[9]= {0.008603, Null}

This would even be faster -- create a matrix using CoeeficientList

In[10]:= Timing[mat = CoefficientList[eq, x];]

Out[10]= {0.002729, Null}

The upshot is that by delaying your computations until the very end, you redo some computations many times over. The lure of using Mathematica functions is that the representation is very compact but for efficiency you sometimes must evaluate the expressions earlier and handle the numbers as matrices. Your profile says you are calling Coefficient too many times so that is a place to start -- create the expression once and use CoeeficientList on the result.

Regards

Neil

Attachments:

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

BTW, if you have n = m =100, then your output will be of order 10^8 elements, is that correct? that needs quite a bit of memory...