Re: Why is Mathematica so slow ?
- To: mathgroup at smc.vnet.net
- Subject: [mg25468] Re: Why is Mathematica so slow ?
- From: "Clifford J. Nelson" <cnelson9 at gte.net>
- Date: Mon, 2 Oct 2000 22:26:50 -0400 (EDT)
- References: <8r1a9q$i50@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Madhusudan Singh wrote:
>
> 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.
a[1] = 3.0
a[2] = 4.0
Creates two variables a[1] and a[2].
a[[1]] = 3.0
a[[2]] = 4.0
assigns 3.0 and 4.0 to the same variable a with different indices, and I
think it is fast in Version 4 and can be compiled with the Compile statement.
Cliff Nelson
>
> (* 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"];