sqrt(64000.0);(32000.00)^(1/4);A straightforward implementation of composite Simpson's rule# composite Simpson's rule
# (straightforward implementation)
compsru := proc( f :: algebraic,
x :: name,
a :: numeric,
b :: numeric,
n :: integer
)
local d, t, s, k;
d := evalf((b-a)/n); # length of a subinterval
t := Vector(n+1);
for k from 1 to n+1 do
t[k] := evalf(a + (k-1)*d)
end do;
s := 0.0;
for k from 1 to n do
s := s + simprule(f,x,t[k],t[k+1])
end do;
return s;
end proc;
#
# Simpson's rule on one interval
#
simprule := proc( f :: algebraic,
x :: name,
a :: numeric,
b :: numeric )
local h, y0, y1, y2;
h := 0.5*(b-a);
y0 := evalf(subs(x=a,f));
y1 := evalf(subs(x=a+h,f));
y2 := evalf(subs(x=b,f));
return (h/3)*(y0 + 4*y1 + y2);
end procread("f:/304/compsru.txt");Testf := sin(x)/x;compsru(f,x,1,5,14);int(f,x=1.0..5.00);solve({c[1]+c[2]=2,c[1]*x[1]+c[2]*x[2]=0,c[1]*x[1]^2+c[2]*x[2]^2=2/3, c[1]*x[1]^3+c[2]*x[2]^3=0},{c[1],c[2],x[1],x[2]});