Quick guide to Mathematica

To get help, place the cursor on the function and press F1
To run all commands in one cell, place the cursor there and press Shift + Enter or Enter on the numerical keyboard.
Output of many commands is suppressed by ; at the end of the command. If you want to see the output, delete ;
Argument of functions must be in square brackets [...], braces {...} are for arrays, ranges etc., use double brackets [[...]] to access elements of arrays
expr /. {x -> a} means substitute a for x in expr, expr // function means apply function to expr, i.e. it is equivalent to function[expr].
D[f[x],x] = derivative of f[x] with respect to x.
To define a function of x one can use f[x_] := 1+(x^2), notice the underscore and the colon, or # nad & like in f := 1+((#)^2) &
N[expr] evaluates expr with machine precision, if you use the decimal point in the expression it will be also evaluated with machine precision
Special characters can be inserted by pressing Esc ... Esc, e.g. Esc p Esc gives π
Clear[...] is used to unset any variables which could have been assigned previously.
Some other useful commands: Simplify, Expand, Factor; Integrate, Series, Sum

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Numerical integration

Newton-Cotes quadrature coefficients

By integrating interpolating polynomials

np=4;
grid = x0+h Range[0,np];
fgrid=Map[f, grid];
points=Table[{grid[[i]],fgrid[[i]]},{i,1,np+1}];
P[x_]=InterpolatingPolynomial[points,x]//Simplify;
Integrate[P[x],{x,x0,x0+np h}]//Simplify//Expand
ω[x_]=Product[x-grid[[i]],{i,1,np+1}];
error=If[Mod[np,2]==0,
Integrate[x ω[x],{x,x0,x0+np h}]/(np+2)!,
Integrate[ω[x],{x,x0,x0+np h}]/(np+1)!
]

By the built-in functions

Needs["NumericalDifferentialEquationAnalysis`"]
NewtonCotesWeights[np+1,0,np]
NewtonCotesError[np+1,f,0,np]

Romberg algorithm (run Basic.Functions.for.TMF057.m first)

Non-periodic function

k=8; (* in the last row, (2^k)+1 points is used *)
a=-1;b=1;
(* f[x_]=Sqrt[1-x^2]; *)
f[x_]=Exp[-x^2];
int=MyRombergIntegrate[f,N[a],N[b],k];
exact=Integrate[f[x],{x,a,b}]
err=int;
Do[Do[If[err[[i,j]]≠0.0,err[[i,j]]=Abs[int[[i,j]]-exact]],{j,1,k+1}],{i,1,k+1}];
MatrixForm[err,15]
Print["Accuracy of the trapezoidal rule with ",(2^k)+1," points is ",N[exact]-int[[k+1,1]]]
Print["Accuracy of the Romberg algorithm with ",(2^k)+1," points is ",N[exact]-int[[k+1,k+1]]]

hg=Table[(b-a)/(2.0^i),{i,0,k}];
data=Table[Table[{hg[[i]],Abs[exact-int[[i,j]]]},{i,j,k+1}],{j,1,5}];
datah=Table[Table[{hg[[i]],((hg[[i]])^j)/6},{i,1,k+1}],{j,{2,4,6,8,10}}];
ListLogLogPlot[{Delete[data,0],Delete[datah,0]},PlotStyle->{Red,Green,Blue,Brown,Orange,Black,Black,Black,Black,Black,Black},Joined->{True},PlotRange->{(10^(-17)),1}]

Periodic function - simple trapezoidal rule is better than Romberg

k=8;
 a=0;b=1;
f[x_]=1/(1+Sin[2π x]/2);
int=MyRombergIntegrate[f,N[a],N[b],k];
(* ScientificForm[int,15] *)
exact=Integrate[f[x],{x,a,b}];
err=int;
Do[Do[If[err[[i,j]]≠0.0,err[[i,j]]=Abs[int[[i,j]]-exact]],{j,1,k+1}],{i,1,k+1}];
MatrixForm[err,15]
Print["Accuracy of the trapezoidal rule with ",(2^k)+1," points is ",N[exact]-int[[k+1,1]]]
Print["Accuracy of the Romberg algorithm with ",(2^k)+1," points is ",N[exact]-int[[k+1,k+1]]]
Plot[f[x],{x,a,b}]

hg=Table[(b-a)/(2.0^i),{i,0,k}];
data=Table[Table[{hg[[i]],Abs[exact-int[[i,j]]]},{i,j,k+1}],{j,1,5}];
datah=Table[Table[{hg[[i]],((hg[[i]])^j)/6},{i,1,k+1}],{j,{2,4,6,8,10}}];
ListLogLogPlot[{Delete[data,0],Delete[datah,0]},PlotStyle->{Red,Green,Blue,Brown,Orange,Black,Black,Black,Black,Black,Black},Joined->{True},PlotRange->{(10^(-16)),1}]

Romberg quadrature coefficients (run Basic.Functions.for.TMF057.m first)

Using Romberg algorithm we can get Simpson’s rule and Boole’s rule (number of points corresponds to behavior of the error) but higher-order quadratures derived from Romberg algorithm have more points than required to obtain given accuracy and thus we get different quadrature rules (less accurate than Newton-Cotes quadrature)

k=2;
int=MyRombergIntegrate[g,α,β,k];
int[[k+1,k+1]]
int[[k+1,k+1]]/.{β->α+(2^k)h}//Simplify

Gauss quadrature

Exaxmple of orthogonal polynomials: Legendre polynomials

LegendreP[3,x]

Plot[Table[LegendreP[n,x],{n,0,5}],{x,-1,1}]

Roots of Legendre polynomials

ListPlot[Flatten[Table[{n,x}/.NSolve[LegendreP[n,x]==0,x],{n,1,40}],1],Axes->False,Frame->True]

General orthogonal polynomials and Gaussian quadrature

Clear[x];
a=-1; b=1;
v[x_]:=1
(* v[x_]:=1/Sqrt[1-x^2] *)
n=3;
MyOrthogonalPolynomials[n,a,b,v, x]
MyGaussianPointsAndWeights[n, a, b, v]

Needs["NumericalDifferentialEquationAnalysis`"]
GaussianQuadratureWeights[n,a,b,16]

Integration with Gaussian quadrature

n=5;
a=-1;b=1;
f[x_]=Sqrt[1-x^2]*Cos[x]; 
v[x_]:=1;(* Gauss-Legendre *)
exact=Integrate[f[x],{x,a,b}]
{xG,wG}=MyGaussianPointsAndWeights[n, a, b, v];
intGL=Sum[wG[[i]]f[xG[[i]]],{i,1,n}]
ScientificForm[intGL-exact,16]
f[x_]=Cos[x];
v[x_]:=Sqrt[1-x^2]; (* Gauss-Chebyshev *)
{xG,wG}=MyGaussianPointsAndWeights[n, a, b, v];
intGC=Sum[wG[[i]]f[xG[[i]]],{i,1,n}]
ScientificForm[intGC-exact,16]