Hello Neil, Thank you for your reply. I agree with you, there is a bug in the calculation and evaluation of BSplineBasis derivatives. Actually, it is correct that the derivatives of BSplineBasis functions involve denominators which are 0. However, this indeterminate result 1/0 is normally replaced by 0. It is a pity that Mathematica does not do it when one uses the following syntax

D[BSplineBasis[{p, U}, i - 1, x], {x, 1}] /. {i -> 1, x -> 0.2}

The workaround you suggest is fine, but it still requires the definitions of derivatives as functions implying definig other functions calling them.... and this slowdown the code. I reformulated the entire code with the overloading operator and in the end the new version is approximately 10 times faster than the initial one. However, it is still much slower than the equivalent matlab version and I am afraid I cannot use it for further developments.

I think I have tried all I could. I don't think that I have other chances I am afraid I will return to Matlab. I will possibly use Mathematica only for small n and m to make comparisons and facilitate the debugging process.

Thank you for all your help!

Enzo

Thank you Neil, I think now I am very close to solve the problem.

Yes, as you siuggested, I do not define any function. I always use = instead of :=. Only in the end, when I need to iterate, I evaluate the symbolic expressions with numbers.

There is only one (I hope minor!) thing that is still not working... very basic issue I guess: can you tell me why I get the error 1/0 indeterminate when I try to evaluate the first derivative of a BsplineBasis? Please see the following short code:

(* Input data *)
p = 2;
n = 5;
r = n + p + 1;
rint = r + 1 - 2 (p + 1);
Uint = Range[0, 1, 1/(rint + 1)]; 
U = N[Join[ConstantArray[0, p], Uint , ConstantArray[1, p]]];

what I do not understand is why the BSplineBasis can be evaluated easily in the following way:

BSplineBasis[{p, U}, i - 1, x] /. {i -> 1, x -> 0.2}

whereas its first derivative CANNOT:

D[BSplineBasis[{p, U}, i - 1, x], {x, 1}] /. {i -> 1, x -> 0.2}

In the example I have chosen to compute the first derivative of the 0th Bspline basis and evaluate it at x=0.2. If you try you'll get "Infinite expression 1/0. encountered. But, as shown in the plot (see file attached), this derivative must DOES exist!

Plot[Evaluate[D[BSplineBasis[{p, U}, 0, x], {x, 1}]], {x, 0, 1}]

Could you please help me on this? After that I think I have done! Many thanks! Enzo

Attachments:

Enzo,

At first glance, I think the Code that Sander posted shows methods to address the issue that I described. Note that he cleverly defines:

bspU[i_]:=bspU[i]=BSplineBasis[{p,U},i-1,u]

This is relying on operator overloading -- for example, when you call the function bspU[3] the calculation is done once so future calls of bspU[3] just return that answer instantly and do no further calculations. If you do this in my simple (contrived) posted example:

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

The first call to define mat would take the full 2 seconds but a subsequent call would take less than a millisecond.

Referring to Sander's code:

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

You should also see if you can define enough of your problem to compute the derivatives once (or do the same function overloading trick). As you iterate over i, the code takes the derivative of bspU[i] repeatedly.

Lastly, Did you try the built in functions for numerically solving PDE's? NSolve and NSolveValue should be helpful since you can express the problem symbolically in Mathematica and have it generate the numerical solution.

Regards

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:

Hello Sander,

thank you very much for your kind reply. I will try all your suggestions (/. intead of //. , SubDivide, etc). However, to answer your main question about “what I want to calculate”, all I need to calculate is reported in the sample code that I posted. My final goal is to build the matrix K (see last line of the code). The only difference in real applications of the full code is that I need to compute many matrices K (say 25) each of them is very big (not like in the posted example). For example, if you change the parameters p, q, n and m by setting them as follows: p = q = 8 and n = m = 100 you will see the code takes ages. Making it UNUSABLE :(

I have only 6 variables which are integers: i, j, ic, jc, r, s. All the rest are real numbers. For me it would be great to treat all variables as reals. How can I tell this to Mathematica? My final result (K) must be made of reals.

Finally, I run the code in debug mode with the command // RuntimeToolsProfile at the end. See profile report in attachment. In this case I set p = q = 2 and n = m = 20. I see that it takes a lot of time in making the substitution Eq //. sub[ic, jc, i, j]. Moreover, I do not understand why it keeps calling many times functions like Coefficient[EqDisc[ic, jc, i, j], {W1[i, j]}] which should be called once for all, instead.

Many thanks,

Enzo

Attachments:

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

Hello Sander,

I went through your reformulation of my example code. It is faster, however I see that you have defined eq = D1w11V2[ij] + D2w112V2[ij]. This works fine but only in this specifc case. I wanted Mathematica to define it because in the general applications of the method the equation can be very long and complex. That's why I moved to Mathematica... to avoid doing it by hand. I hope that after my explanation of the whole meaning of the code (see my yestrday post) you have a clearer picture of my final goal.

Many thaks,

Enzo