--  Romberg integration
--  This is obtained by applying extrapolation to the trapezoid rule.
--  by Michael E. Gage      gage@cc.rochester.edu

-- ( use Radians! )

--Define function
f(x)=ct:=ct+1,sin(x)
start=0
end=X

--Definition of Romberg extrapolation scheme
Romberg(j,k) = trapezoid(2^(j-1)) when k=1,
               (4^(k-1)*Romberg(j,k-1)-
               Romberg(j-1,k-1))/(4^(k-1)-1) when 1<k and k<=j

--Definition of trapezoid rule
trapezoid(n)=delt:=(end-start)/n,
delt*(sum(f(start+k*delt),k,1,n-1)+0.5*f(start)+0.5*f(end));

-- plot the error
Xmin=0; Xmax=2*pi
ct:=0:
plot Romberg(4,3) - (-cos(X)+1)
ct:2222.000000

ct:=0:
plot trapezoid(20) - (-cos(X)+1)
ct:2121.000000