%高斯系列公式数值积分
%格式为:q=IntGauss(f,a,b,n,AK,XK)
function q=IntGauss(f,a,b,n,AK,XK)
% 被积函数:f
%积分区间:(a,b)
%所采用的高斯积分点个数:n;
%自定义系数:AK;
%自定义积分点:XK;
%积分值:q
if (n<5 && nargin==4)
AK=0;
XK=0;
else %如果n>4,则节点和系数由调用者给出
XK1=((b-a)/2*XK+(a+b)/2);
q=((b-a)/2)*sum(AK.*subs(sym(f),findsym(f),XK1));
end
ta = (b-a)/2;
tb = (a+b)/2;
switch n
case 0,
q=2*ta*subs(sym(f),findsym(sym(f)),tb);
case 1,
q=ta*(subs(sym(f),findsym(sym(f)),ta*0.5773503+tb)+...
subs(sym(f),findsym(sym(f)),-ta*0.5773503+tb));
case 2,
q=ta*(0.55555556*subs(sym(f),findsym(sym(f)),ta*0.7745967+tb)+...
0.55555556*subs(sym(f),findsym(sym(f)),-ta*0.7745967+tb)+...
0.88888889*subs(sym(f),findsym(sym(f)),tb));
case 3,
qq=ta*(0.3478548*subs(sym(f),findsym(sym(f)),ta*0.8611363+tb)+...
0.3478548*subs(sym(f),findsym(sym(f)),-ta*0.8611363+tb)+...
0.6521452*subs(sym(f),findsym(sym(f)),ta*0.3398810+tb)...
+0.6521452*subs(sym(f),findsym(sym(f)),-ta*0.3398810+tb));
case 4,
q=ta*(0.2369269*subs(sym(f),findsym(sym(f)),ta*0.9061793+tb)+...
0.2369269*subs(sym(f),findsym(sym(f)),-ta*0.9061793+tb)+...
0.4786287*subs(sym(f),findsym(sym(f)),ta*0.5384693+tb)...
+0.4786287*subs(sym(f),findsym(sym(f)),-ta*0.5384693+tb)+...
0.5688889*subs(sym(f),findsym(sym(f)),tb));
end