[Date Index]
[Thread Index]
[Author Index]
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
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.
>
> (* 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**
| |