Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2000
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

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




  • Prev by Date: Point inside a plygon (again)
  • Previous by thread: Point inside a plygon (again)