matlab quad2,matlab中quad.m函数,我把它简化了,function [Q] = quad2(funfcn,a,b)f = fcnchk(funfcn);h = 0.13579*(...

matlab中quad.m函数,我把它简化了,function [Q] = quad2(funfcn,a,b)f = fcnchk(funfcn);h = 0.13579*(b-a);x = [a a+h a+2*h (a+b)/2 b-2*h b-h b];y = f(x);[Q(1)] = ...quadstep(f,x(1),x(3),y(1),y(2),y(3));[Q(2)] = ...quadstep(f,x(3),x(5),y(3),y(4),y

2016-03-30

matlab中quad.m函数,我把它简化了,

function [Q] = quad2(funfcn,a,b)

f = fcnchk(funfcn);

h = 0.13579*(b-a);

x = [a a+h a+2*h (a+b)/2 b-2*h b-h b];

y = f(x);

[Q(1)] = ...

quadstep(f,x(1),x(3),y(1),y(2),y(3));

[Q(2)] = ...

quadstep(f,x(3),x(5),y(3),y(4),y(5));

[Q(3)] = ...

quadstep(f,x(5),x(7),y(5),y(6),y(7));

Q = sum(Q);

function [Q] = quadstep(f,a,b,fa,fc,fb)

tol = 1.e-6;

h = b - a;

c = (a + b)/2;

x = [(a + c)/2,(c + b)/2];

y = f(x);

fd = y(1);

fe = y(2);

% Three point Simpson's rule.

Q1 = (h/6)*(fa + 4*fc + fb);

% Five point double Simpson's rule.

Q2 = (h/12)*(fa + 4*fd + 2*fc + 4*fe + fb);

% One step of Romberg extrapolation.

Q = Q2 + (Q2 - Q1)/15;

if abs(Q2 - Q)