Gaussian Quadrature

GaussianQuadWeights[n_,a_,b_]:=With[{ruledata=NIntegrate`GaussRuleData[n,MachinePrecision]},Transpose[{Rescale[ruledata[[1]],{0,1},{a,b}],Abs[b-a]ruledata[[2]]}]];

GQuadPlot[f_,{a_,b_},n_,max_ :20]:=Module[{int,approxint,qp,data,lp,weights,frst,lst,polys,ep,tp},
With[{x=x,i=i},
int = N[Integrate[f[x],{x,a,b}]];
weights = GaussianQuadWeights[n,a,b];
approxint = Total[f[First[#]]*Last[#]&/@weights];
frst  = {First[weights[[1]]],0};
lst = {First[weights[[n]]],0};
polys={Hue[Mod[((GoldenRatio(#[[1,2]] + #[[2,2]]))/2),1]],Polygon[{{#[[1,1]],0},Sequence@@#,{#[[2,1]],0}}]}&/@Partition[weights,2,1];

lp =ListPlot[weights,Filling->{1->{Axis,Black}},PlotStyle->{Black,PointSize[Large]},PlotRange->{{a,b},{0,Max[Last/@weights]}},Prolog->polys,Axes->{True,False},Ticks->{
Part[Transpose[{First/@weights,Subscript[x,#]&/@Range[1,n]}]
,{1,Sequence @@Range[3,n-2,2],n}]
,Automatic},ImageSize->{250,170},
PlotRangePadding->.2 Max[Last/@weights],Epilog->{{Red,PointSize[Large],Point[{First[#],0}]&/@weights},Part[MapIndexed[Text[Subscript[TraditionalForm[w],ToString[First[#2]]],#1,{0,-2}]&,weights],{1,Sequence @@Range[3,n-2,2],n}]}];

qp = Plot[f[x],{x,a,b},Filling->{1->{Axis,Lighter[Green,.2]}},Epilog->{{Black,PointSize[Large],({Line[{{First[#],0},{First[#],f[First[#]]}}],Point[{First[#],f[First[#]]}]}&/@weights)},{Red,PointSize[Large],Point[{First[#],0}]&/@weights}},ImageSize->{250,170}];

data =Table[{i,GQuad[f,{a,b},i]},{i,2,n-1}];
AppendTo[data,{n,approxint}];

ep = ListPlot[data,Epilog->{Line[{{0,int},{max,int}}],Text["true integral",{max,int},{1,-1}]},ImageSize->{250,190},PlotRange->{{1,max},
{int-.3*Abs[int],int + .3*Abs[int]}},AxesOrigin->{1,0},PlotRangePadding->{0.2,.1*Abs[int]},FrameTicks->{Range[0,max,2],Automatic},Frame-> {True,True,False,False},FrameLabel->{"number of points","estimate"},RotateLabel->True,PlotStyle->{Lighter[Green,.2],PointSize[Medium]},Filling->{1->{Bottom,Black}},ImageSize->{250,170}];

tp = 
Pane[Text@Style[Column[{"    ",
Row[{"   ",Underoverscript["∑",Row[{Style["i",Italic],"=","1"}],ToString[n]],Row[{Row[{"(",TraditionalForm[f[(Subscript[x, i])]],")"}],TraditionalForm[(Subscript[w, i])]}]," = ",approxint}],Row[{" ≃ ",Subsuperscript["∫",a,b],"(",TraditionalForm[f[x]],")",Row[{"\[DifferentialD]",Style["x",Italic]}]}],Row[ {" = ",TraditionalForm@N[Integrate[f[x],{x,a,b}],20]}]}],20],330];

Grid[{{
Text@Column[
{Column[{"integrand",qp},Center],Column[{"weights and sampling points",lp},Center]},Center],
Column[{tp,ep},Center,8.3]}}]]

];

GQuad[f_,{a_,b_},n_]:=Total[f[First[#]]*Last[#]&/@GaussianQuadWeights[n,a,b]];

funcandrange = {
{Function[x,3 Cos[x] + Sin[x]],{-2,5/2}},
{Function[x,Exp[x] - x^2],{-2,2}},
{Function[x,x^2 -1],{-2,2}},
{Function[x,Cos[9 x] + 5x ],{-2,2}},
{Function[x,Sin[12 x] - x^2],{-1,2}}
};
functions = Rule[Evaluate[TraditionalForm[Last[First[#]]]],#]&/@ funcandrange;
GetFunc[str_]:= First[str /. functions];
GetRange[str_]:=Last[str /. functions];

Manipulate[GQuadPlot[GetFunc[str] ,GetRange[str],n],{{n,2,"number of points"},2,25,1,Appearance->"Labeled"},{{str,functions[[1,1]],"integrand"},First/@ functions,ControlType->PopupMenu},TrackedSymbols:>{n,str},SaveDefinitions->True]

CAPTION

Gaussian quadrature with n points gives exact integral estimates for polynomials of degree 2 n-1 or less. For smooth functions, increasing the number of points of the Gaussian rule improves the integral estimate.

THIS NOTEBOOK IS THE SOURCE CODE FROM

"Gaussian Quadrature" from the Wolfram Demonstrations Project
 http://demonstrations.wolfram.com/GaussianQuadrature/

Contributed by: Chris Maes

Additional contributions by: Anton Antonov

A full-function Wolfram Mathematica system (Version 6 or higher) is required to edit this notebook.
GET WOLFRAM MATHEMATICA »

© Wolfram Demonstrations Project & Contributors   |   Terms of Use   |   Make a new version of this Demonstration »