Re: Why is Mathematica so slow ?

• To: mathgroup at smc.vnet.net
• Subject: [mg25781] Re: Why is Mathematica so slow ?
• From: Robert Knapp <rknapp at wolfram.com>
• Date: Wed, 25 Oct 2000 03:53:42 -0400 (EDT)
• Organization: Wolfram Research, Inc.
• References: <8r1a9q\$i50@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```If I surmise correctly that the main complexity in your code is solving
a tridiagonal linear system, then it can be greatly speeded up, even if
you cannot use the Mathematica SparseLinearSolve command.

I showed how to do this in a talk last year.  You can find the notebook
at

http://members.wri.com/rknapp/Talks/DeveloperConference99/Tridiagonal.nb

>
> 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,
>
>
> (* 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: Re: Summary: List element manipulation
• Next by Date: Re: How do you create a ListContourPlot from a mesh of non-equally spaced real numbers?
• Previous by thread: Re: Why is Mathematica so slow ?
• Next by thread: NDSolve with a function which calls N