Why is Mathematica so slow ?
- To: mathgroup at smc.vnet.net
- Subject: [mg25415] Why is Mathematica so slow ?
- From: Madhusudan Singh <chhabra at eecs.umich.edu>
- Date: Fri, 29 Sep 2000 01:07:13 -0400 (EDT)
- Organization: EECS Dept. Univ. of Michigan
- Sender: owner-wri-mathgroup at wolfram.com
Hi
I am solving a problem in a numerical linear algebra course that
involves solving a two point boundary value problem with finite element
method and LU factorisation. Since this is a part of a course, I am not
permitted to use inbuilt functions in Mathematica.
The matrix generated is a tridiagonal matrix. It is sparse. The
order varies from 1 to 2^14. I solved this problem by using the
following code (end of posting).
I did the same problem on C (since I am using Mathematica purely as
a simple programming language, it does not matter if I use C instead.)
To solve for matrix orders 1 through 2^14 (14 steps) it takes about 4
minutes to solve the problem with C, and more than 7 hours with
Mathematica. The two codes are functionally identical. What is going on
?
With regards,
Madhusudan Singh.
(* The definition and initialisation section.*)
Clear["'*"]; Off[
Part::"pspec"]; Off[General::"spell1"]; Off[Part::"partw"];
a0 := 0;
b0 := 1;
pcons := 14;
g[x_] := x^2;
f[x_] := (1 + 4 x + 2 x^2 - x^4) Exp[x];
phi[x_] := (1 - x^2) Exp[x];
h[p_] := 1/(2^(p));
n[p_] := (b0 - a0)/h[p] - 1;
Array[maxnorm, pcons];
Print["h(p) ||uh-phih|| ||uh-phih||/h2"];
(*The loop over all p's."*)
Do[{
Clear[meshsize]; Clear[number];
meshsize = h[p];
number = n[p];
Clear[a]; Clear[b]; Clear[c];
Array[a, {number}];
Array[b, {number}];
Array[c, {number}];
Clear[diffupper]; Clear[difflower]; Clear[diffdiag];
(*Definition of the three diagonals.
The matrix definition is eschewed as the matrix is sparse
and \
(for larger p's) memory can become an issue.*)
Do[{a[k] := N[2/(meshsize)^(2) + g[a0 + k meshsize]];
b[k + 1] := N[-1/(meshsize)^(2)];
c[k] := N[-1/(meshsize)^(2)];}, {k, 1, number - 1}];
a[number] := N[2/(meshsize)^(2) + g[a0 + number meshsize]];
diffdiag := Table[a[k], {k, 1, number}];
diffupper := Table[c[k], {k, 1, number - 1}];
difflower := diffupper;
Clear[alpha]; Clear[beta]; Clear[gamma];
Array[gamma, number - 1];
Array[alpha, number];
Array[beta, number]; i := 1;
(*Calculation of alpha,
beta and gamma parameters from the original eigenvalue
equation.*)
\
Do[{gamma[i] := diffupper[[i]]; }, {i, 1, number - 1}];
beta[1] = 0;
beta[2] = difflower[[1]]/diffdiag[[1]];
alpha[1] = diffdiag[[1]];
Do[{alpha[i] = diffdiag[[i]] - diffupper[[i - 1]] beta[i];
beta[i + 1] = difflower[[i]]/alpha[i];}, {i, 2, number - 1}];
alpha[number] =
diffdiag[[number]] - diffupper[[number - 1]] beta[number];
(*Forward elimination.*)
Clear[y]; Clear[x];
Array[y, number];
y[1] = f[a0 + 1 meshsize] + 1/(meshsize)^(2);
Do[{y[i] = f[a0 + i meshsize] - beta[i] y[i - 1];}, {i, 2,
number}];
(*Backward substitution.*)
Clear[ulist];
Array[x, number];
x[number] = y[number]/alpha[number];
Do[{x[i] = (y[i] - gamma[i] x[i + 1])/alpha[i];}, {i, number - 1,
1, -1}];
ulist := Table[Abs[x[i] - phi[a0 + i meshsize]], {i, 1, number}];
maxnorm[p] = Max[ulist];
resultstream =
OpenAppend["result13.dat", FormatType -> OutputForm,
PageWidth -> Infinity];
Write[resultstream, N[meshsize], " ", maxnorm[p], " ",
maxnorm[p]/(meshsize)^(2)];
Print[N[meshsize], " ", maxnorm[p], " ",
maxnorm[p]/(meshsize)^(2)];
Close[resultstream];}, {p, 1, pcons}];
(*Plotting the results*)
maxnormlist := Table[{h[p], maxnorm[p]}, {p, 1, pcons}];
algorithmefficiency := Table[{h[p], maxnorm[p]/(h[p])^(2)}, {p, 1,
pcons}];
Display["maxnorm.eps", ListPlot[maxnormlist, PlotJoined -> True],
"EPS"];
Display["algoeff.eps", ListPlot[algorithmefficiency, PlotJoined ->
True],
"EPS"];