gauss_legendre2.m

% 任意个数插值节点的高斯-勒让德求积公式:
% 思路: 其内涵还是n阶拉格朗日插值! 只不过原先的拉格朗日插值点都是等间距均匀分布的那种,
%       高斯-勒让德就是点数一定时, 挑选那些能让精度提高的点(选最好的点)!而不再使用无特殊意义的等间距分隔点。
%       补充: 高斯-勒让德"不能"与前面的方法混用(比如选最好点后再龙贝格)!
%             因为前面所有的方法要求"必须是等间距"的点才能写出"那样的公式"!
% 优点: 小区间/积分限上,限制了只能有几个控制点时, 用高斯积分效果不错! 
% 缺点: 大区间上, 高斯-勒让德的求积节点必须得增多! 因为内涵还是n阶拉格朗日插值!——这个缺点非常大!
% 建议: 没有控制点限制时, 我还是推荐用龙贝格! 因为它是万能适用的高精度! 

clear; clc;
format long;

syms x;
n = double(input('输入使用几点(n)的高斯-勒让德插值:'));
% n点插值的高斯-勒让德多项式P:
f = x^2 -1;
fprintf('%d点高斯-勒让德多项式为:\n',n)
P = vpa(1/(2^n*factorial(n)) * diff(f^n,x,n))
% n点高斯-勒让德插值节点Xi:
Xi = sort(double(solve(P)))';

% n点高斯-勒让德插值节点对应的插值系数Ai:
xnum = Xi;
% l用来记录拉格朗日多项式的: n点拉格朗日插值多项式有n个基函数! 
% 高斯-勒让德求积系数就是对"每个多项式/基函数"求其[-1 1]定积分——
% 插值点数、多项式/基函数个数、求积系数个数一样，相互对应。
l = sym(zeros(1,n));  
Ai = zeros(1,n);
for m = 1:n
    l(m) = prod(x - xnum([1:m-1 m+1:n]))/prod(xnum(m) - xnum([1:m-1 m+1:n])); 
    Ai(m) = double(int(l(m),x,-1,1));
end

fprintf('%d点高斯-勒让德插值节点为:\n',n);
Xi
fprintf('%d点高斯-勒让德插值节点对应的系数为:\n',n);
Ai


% 对具体函数的积分近似: 
syms t;
up = double(input('输入积分上限:'));
low = double(input('输入积分下限:'));
% 积分限x必须是[-1 1], 因此这里统一做一下转换, 即任何限都可以:
t = (up-low)/2*x + (up+low)/2;
y = 4/(1+t^2);  % 每次修改这里的函数即可, 注意这里自变量是t (修改√√√)

Result = 0;  % 记录最终近似结果
for num = 1:n
    x = Xi(num);
    y_tmp = (up-low)/2*double(subs(y));
    Result = Result + Ai(num)*y_tmp;
end

syms w;
y_r = 4/(1+w^2);  % 实际的结果, 注意这里自变量是w  (修改√√√)
R = int(y_r,w,low,up);
fprintf('%d点高斯-勒让德积分近似结果为:%.10f\n', n, Result);
fprintf('真实结果为:%.10f\n', R);