MATLAB实现区间[a,b]上的Gauss-Legendre积分

Question

区间[a,b]上2点Gauss-Legendre积分如下，计算[0,π/4]上cosx的积分结果和精确值相差太大，如何修改？

%区间[a,b]上2点Gauss-Legendre积分
function y=GLdi2ab(fun,a,b)
f1=fun((b+a)/2+((b-a)/2)*(-1/sqrt(3)));
f2=fun((b+a)/2+((b-a)/2)*(1/sqrt(3)));
y=((b-a)*2)*(f1+f2);
fprintf("该区间上2点Gauss-Legendre积分积分公式的结果为\n");
disp(y);

区间[a,b]上3点Gauss-Legendre积分如下，计算[0,π/4]上cosx的积分结果和精确值相差太大，如何修改？

function y=GLdi3ab(fun,a,b)
f1=fun((b+a)/2+((b-a)/2)*(-sqrt(15)/5));
f2=fun((b+a)/2);
f3=fun((b+a)/2+((b-a)/2)*(sqrt(15)/5));
y=((b-a)*2)*((5/9)*f1+(8/9)*f2+(5/9)*f3);
fprintf("该区间上3点Gauss-Legendre积分积分公式的结果为\n");
disp(y);

计算误差的函数如下

function e=merror(fun1,fun2,a,b,fun)
%fun1为数值积分计算公式
%fun2为y=cosx的定积分计算
%fun为数值积分的被积函数
%a,b分别为积分区间的左端点，右端点
y1=fun1(@fun,a,b);
y2=fun2(a,b);
e=abs(y1-y2);
fprintf("误差为\n");
disp(e)

Answer 1

🕗晚上好，题主！目前：🟢[在线]
📢参考通义千问、由郭老师整理提供！
📣今日诗号：古岂无人，意行千里论春秋；山月笑卧，孤标尘外一剑宿。

首先，你的GLdi2ab和GLdi3ab函数中使用的是Gauss-Legendre积分公式，但是它们并不是在[a, b]上精确地应用。对于2点和3点的Gauss-Legendre积分，我们需要找到相应的节点（也称为Gauss点）和权重。这些节点和权重是根据Gauss-Legendre多项式来确定的，而不是简单地将积分区间等分为两段或三段。

对于2点Gauss-Legendre积分，节点和权重分别为：

• 节点：( x_1 = \frac{1}{2} - \frac{\sqrt{3}}{6} ) 和 ( x_2 = \frac{1}{2} + \frac{\sqrt{3}}{6} )
• 权重：( w_1 = \frac{1}{2} ) 和 ( w_2 = \frac{1}{2} )

对于3点Gauss-Legendre积分，节点和权重分别为：

• 节点：( x_1 = \frac{-\sqrt{15}}{5} )，( x_2 = 0 )，( x_3 = \frac{\sqrt{15}}{5} )
• 权重：( w_1 = \frac{5}{9} )，( w_2 = \frac{8}{9} )，( w_3 = \frac{5}{9} )

修改后的代码如下：

function y=GLdi2ab(fun,a,b)
    x1 = (1/2) - (sqrt(3)/6);
    x2 = (1/2) + (sqrt(3)/6);
    w1 = 1/2;
    w2 = 1/2;

    f1 = fun(x1*(b-a)+a);
    f2 = fun(x2*(b-a)+a);
    y = (b-a) * (w1*f1 + w2*f2);
    fprintf("该区间上2点Gauss-Legendre积分积分公式的结果为\n");
    disp(y);
end

function y=GLdi3ab(fun,a,b)
    x1 = (-sqrt(15))/5;
    x2 = 0;
    x3 = (sqrt(15))/5;
    w1 = 5/9;
    w2 = 8/9;
    w3 = 5/9;

    f1 = fun(x1*(b-a)+a);
    f2 = fun(x2*(b-a)+a);
    f3 = fun(x3*(b-a)+a);
    y = (b-a) * (w1*f1 + w2*f2 + w3*f3);
    fprintf("该区间上3点Gauss-Legendre积分积分公式的结果为\n");
    disp(y);
end

至于merror函数，它已经正确计算了两个函数的积分结果之间的误差。如果fun1和fun2分别代表数值积分计算公式和精确积分计算，那么这个函数没有问题。

为了计算cos(x)在[0, π/4]上的积分，你需要提供一个fun函数，例如：

function y = cos_x(x)
    y = cos(x);
end

然后调用merror函数：

fun1 = @GLdi2ab;
fun2 = @(x) integrate(@(t) cos(t), 0, pi/4); % 使用integral函数计算精确值
a = 0;
b = pi/4;
e = merror(fun1, fun2, a, b, @cos_x);
fprintf("误差为\n");
disp(e);

这样应该能得到更接近精确值的积分结果。如果你仍然发现误差很大，请检查fun函数是否正确实现了cos(x)。