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

Re: Effect of Simplify on numeric vs symbolic

• To: mathgroup at smc.vnet.net
• Subject: [mg56027] Re: Effect of Simplify on numeric vs symbolic
• From: Paul Abbott <paul at physics.uwa.edu.au>
• Date: Thu, 14 Apr 2005 08:54:18 -0400 (EDT)
• Organization: The University of Western Australia
• References: <d3g8vq$t01$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

In article <d3g8vq$t01$1 at smc.vnet.net>, carlos at colorado.edu wrote:

> In my previous example of FEM symbolic vs numeric speed all
> computations were with numbers.  No variables were involved.
> The difference (minutes vs seconds) due largely due to Simplify.

I could find no mention of Gauss-Legendre quadrature in the index of 
your INTRODUCTION to FINITE ELEMENT METHODS course at 

  http://caswww.colorado.edu/courses.d/IFEM.d/Home.html ?
 
> Here is a simpler example that illustrates the point, also taken from
> the same course.  LineGaussRuleInfo returns weights and abcissae of 1D
> Gauss-Legendre quadrature rules of 1 to 5 points. Argument numer asks
> for exact information if False; if True it returns N[information].
> All computations involve nunbers only.

The expressions involve arithmetic with complicated radicals (that can 
be avoided -- see below). Manipulation of such expression is often very 
slow.
 
> Results of calls (at end of post) on a Mac laptop:
> 
> numer rule   n  Eval time  Simplify time
> True    4   50   0.00 sec     0.00 sec
> False   4   50   0.00 sec    78.01 sec  (* Note: FullSimplify needed *)
> False   4   50   0.00 sec     0.00 sec  (* found result in cache *)

I ran your code (Mac OS X 10.3.8, Mathematica 5.1) and don't get these 
timings? The _first_ time I call

 PerverseExpression[4, 50, False]

I get the correct answer in 0.16 sec.

> LineGaussRuleInfo[{rule_,numer_},point_]:= Module[
>   {g2={-1,1}/Sqrt[3],w3={5/9,8/9,5/9},
>    g3={-Sqrt[3/5],0,Sqrt[3/5]},
>    w4={(1/2)-Sqrt[5/6]/6, (1/2)+Sqrt[5/6]/6,
>        (1/2)+Sqrt[5/6]/6, (1/2)-Sqrt[5/6]/6},
>    g4={-Sqrt[(3+2*Sqrt[6/5])/7],-Sqrt[(3-2*Sqrt[6/5])/7],
>         Sqrt[(3-2*Sqrt[6/5])/7], Sqrt[(3+2*Sqrt[6/5])/7]},
>    g5={-Sqrt[5+2*Sqrt[10/7]],-Sqrt[5-2*Sqrt[10/7]],0,
>         Sqrt[5-2*Sqrt[10/7]], Sqrt[5+2*Sqrt[10/7]]}/3,
>    w5={322-13*Sqrt[70],322+13*Sqrt[70],512,
>        322+13*Sqrt[70],322-13*Sqrt[70]}/900,
>    i=point,p=rule,info={{Null,Null},0}},
>   If [p==1, info={0,2}];
>   If [p==2, info={g2[[i]],1}];
>   If [p==3, info={g3[[i]],w3[[i]]}];
>   If [p==4, info={g4[[i]],w4[[i]]}];
>   If [p==5, info={g5[[i]],w5[[i]]}];
>   If [numer, Return[N[info]], Return[info] ];
> ];

The abscissas and weights for arbitrary n can be computed directly as 
follows:

 x[n_, m_] := Root[Function[x, LegendreP[n, x]], m]

 w[n_, m_] := 2 (1 - x[n, m]^2)/((n + 1)^2 LegendreP[n + 1, x[n, m]]^2)

a direct implementation of equation (10) at

 http://mathworld.wolfram.com/Legendre-GaussQuadrature.html

To convert these Root objects to radicals, we use ToRadicals:

 TableForm[ToRadicals[Table[x[n, m], {n, 2, 5}, {m, 1, n}]]]

Similarly, the weights reduce to

 TableForm[RootReduce[Table[w[n, m], {n, 3, 5}, {m, 1, n}]]]

These results agree, as expected with your gn's and wn's.

However, there really is no need to convert the abscissas and weights 
expressions to radicals. In general, especially for large n (in fact, 
even for n=6), it is better to work with representations involving Root 
objects and use RootReduce to simplify combinations of such objects.

> PerverseExpression[rule_,n_,numer_]:=Module[{xw,tEval,tSymb,x},
>   xw=Table[LineGaussRuleInfo[{rule,numer},i],{i,1,rule}];
>   {tEval,x}=Timing[Product[(i+2)*xw[[i,1]]^n,{i,1,rule}]];
>   {tSymp,x}=Timing[FullSimplify[x]];
>   Return[{tEval,tSymp,x}]];

This expression is just 

  g[n_, p_] := Product[(i + 2) x[n, i]^p, {i, 1, n}]

Computing 

  RootReduce[g[4, 50]]

takes about 0.4 sec and gives the same result as

  Print[PerverseExpression[4,50,False]]

Cheers,
Paul

-- 
Paul Abbott                                   Phone: +61 8 6488 2734
School of Physics, M013                         Fax: +61 8 6488 1014
The University of Western Australia      (CRICOS Provider No 00126G)         
35 Stirling Highway
Crawley WA 6009                      mailto:paul at physics.uwa.edu.au 
AUSTRALIA                            http://physics.uwa.edu.au/~paul