复化梯形公式，Newton-Cotes公式，变量代换后的复化梯形公式，Gauss-Legendre公式，Gauss-Jacobi公式插值积分的精确度比较

程序员文章站 2022-07-12 14:15:14

...

1.问题

分别计算积分
I c = ∫ 0 1 cos ⁡ x x d x = 1.809048475800... I_c=\int_0^1\frac{\cos{x}}{\sqrt{x}}dx=1.809048475800... Ic=∫01x cosxdx=1.809048475800...

I s = ∫ 0 1 sin ⁡ x x d x = 0.620536603446 I_s=\int_0^1\frac{\sin{x}}{\sqrt{x}}dx=0.620536603446 Is=∫01x sinxdx=0.620536603446

根据和真实值的对比，探讨各种方法的优劣。

2.数学理论和算法

1. n个等分区间的复化梯形公式n=100 : 100 : 1000

T n = ∑ k = 0 n − 1 ( x k + 1 − x k ) 2 ( f ( x k ) + f ( x k + 1 ) ) T_n=\sum_{k=0}^{n-1}\frac{(x_{k+1}-x_k)}{2}(f(x_k)+f(x_{k+1})) Tn=k=0∑n−12(xk+1−xk)(f(xk)+f(xk+1))

2. 带权的Newton-Cotes公式n=100 : 100 : 1000

推导过程：
f ( x ) = 1 , i n t e g r a l = 2 = A + B ; f ( x ) = x , i n t e g r a l = 2 / 3 = B ; A = 4 / 3 , B = 2 / 3 f(x)=1,integral=2=A+B; f(x)=x,integral=2/3=B; A=4/3,B=2/3 f(x)=1,integral=2=A+B;f(x)=x,integral=2/3=B;A=4/3,B=2/3

思路：对区间[0,h],由公式(3)代数精度为1,可得A=(4/3)(h^{(1/2));B=(2/3)(h}(1/2))
x i = a + i h , i = 0 , 1 , ⋅ ⋅ ⋅ ， n , h = b − a n x_i=a+ih,i=0,1,\cdot\cdot\cdot，n,h=\frac{b-a}{n} xi=a+ih,i=0,1,⋅⋅⋅，n,h=nb−a

P n ( x ) = ∑ i = 0 n ω ( x ) ( x − x i ) ω ′ ( x ) f ( x i ) , ω ( x ) = ( x − x 0 ) ( x − x 1 ) ⋅ ⋅ ⋅ ( x − x n ) P_n(x)=\sum_{i=0}^n\frac{\omega(x)}{(x-x_i)\omega'(x)}f(x_i), \omega(x)=(x-x_0)(x-x_1)\cdot\cdot\cdot(x-x_n) Pn(x)=i=0∑n(x−xi)ω′(x)ω(x)f(xi),ω(x)=(x−x0)(x−x1)⋅⋅⋅(x−xn)

∫ a b f ( x ) d x ≈ ∫ a b P n ( x ) d x = ∫ a b ( ∑ i = 0 n ω ( x ) ( x − x i ) ω ′ ( x i ) f ( x i ) ) d x = ∑ i = 0 n ( ∫ a b ω ( x ) ( x − x i ) ω ′ ( x i ) d x ) f ( x i ) = ∑ i = 0 n A i f ( x i ) \int_a^bf(x)dx\approx\int_a^bP_n(x)dx=\int_a^b(\sum_{i=0}^n\frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i))dx=\sum_{i=0}^n(\int_a^b\frac{\omega(x)}{(x-x_i)\omega'(x_i)}dx)f(x_i)=\sum_{i=0}^nA_if(x_i) ∫abf(x)dx≈∫abPn(x)dx=∫ab(i=0∑n(x−xi)ω′(xi)ω(x)f(xi))dx=i=0∑n(∫ab(x−xi)ω′(xi)ω(x)dx)f(xi)=i=0∑nAif(xi)

A i = ∫ a b ω ( x ) ( x − x i ) ω ′ ( x i ) d x A_i=\int_a^b\frac{\omega(x)}{(x-x_i)\omega'(x_i)}dx Ai=∫ab(x−xi)ω′(xi)ω(x)dx

3. 变量代换后的n个等分区间的复化梯形公式n=20 : 20 : 200

x = t 2 , I c ′ = ∫ 0 1 2 cos ⁡ ( t 2 ) d t , I s ′ = ∫ 0 1 2 sin ⁡ ( t 2 ) d t x=t^2,I_c'=\int_0^12\cos(t^2)dt,I_s'=\int_0^12\sin(t^2)dt x=t2,Ic′=∫012cos(t2)dt,Is′=∫012sin(t2)dt

复化梯形公式如1.

4. Gauss-Legendre积分n=1 : 1 : 5

t = 0 + 1 2 + 1 − 0 2 T = 1 + T 2 t=\frac{0+1}{2}+\frac{1-0}{2}T=\frac{1+T}{2} t=20+1+21−0T=21+T

I c = ∫ − 1 1 cos ⁡ ( 1 + T 2 ) 2 d T I_c=\int_{-1}^1\cos(\frac{1+T}{2})^2dT Ic=∫−11cos(21+T)2dT

I s = ∫ − 1 1 sin ⁡ ( 1 + T 2 ) 2 d T I_s=\int_{-1}^1\sin(\frac{1+T}{2})^2dT Is=∫−11sin(21+T)2dT

5. 由正交多项式推出相应Gauss型积分公式n=1 : 1 : 5

3.程序与结果

1. 复化梯形公式

编辑器：

syms x
f(x)=cos(x)/sqrt(x);
g(x)=sin(x)/sqrt(x);
errof=[];
errog=[];
for n=100:100:1000
    h=1/n;
    x=linspace(1/n,1,n);
    Tf=eval((h*(0+sum(f(x))-f(x(1,n))/2)))
    Tg=eval((h*(0+sum(g(x))-g(x(1,n))/2)))
    errof(n/100)=abs(Tf-Tfbar);
    errog(n/100)=abs(Tg-Tgbar);
end
Tfbar=eval(int(f(x),0,1))
Tgbar=eval(int(g(x),0,1))
errof,errog

Tf=1.663004
Tg=0.620330
Tf=1.705784
Tg=0.620463
Tf=1.724734
Tg=0.620497
Tf=1.736030
Tg=0.620511
Tf=1.743739
Tg=0.620518
Tf=1.749429
Tg=0.620522
Tf=1.753852
Tg=0.620525
Tf=1.757417
Tg=0.620527
Tf=1.760370
Tg=0.620530
Tf=1.762870
Tg=0.620530
Tfbar=1.809048
Tgbar=0.620537
errof=
  列 1 至 2 
   0.146045   0.103265
  列 3 至 4
   0.084315   0.073018
  列 5 至 6
   0.065309   0.059619
  列 7 至 8
   0.055196   0.051631
  列 9 至 10
   0.048679   0.046181
errog=
   1.0e-03 *
  列 1 至 2
   0.206890   0.073250
  列 3 至 4
   0.039897   0.025923
  列 5 至 6
   0.018554   0.014117
  列 7 至 8
   0.011204   0.009172
  列 9 至 10
   0.007687   0.006563

2. 带权的Newton-Cotes公式

编辑器：

syms x
f(x)=cos(x)/sqrt(x);
g(x)=sin(x)/sqrt(x);
errof=[];
errog=[];
for n=100:100:1000
    h=1/n;
    x=linspace(1/n,1,n);
    Tf=eval((4/3)*n^(-1/2)+(2/3)*n^(-1/2)*cos(1/n)+(h*(f(x(1,1))/2+sum(f(x))-f(x(1,n))/2)))
    Tg=eval((2/3)*n^(-1/2)*sin(1/n)+(h*(g(x(1,1))/2+sum(g(x))-g(x(1,n))/2)))
    errof(n/100)=abs(int(f(x),0,1))
    errog(n/100)=abs(int(g(x),0,1))
end
Tfbar=eval(int(f(x),0,1))
Tgbar=eval(int(g(x),0,1))
errof,errog

Tf=1.912998
Tg=0.621496
Tf=1.882559
Tg=0.620875
Tf=1.869071
Tg=0.620721
Tf=1.861030
Tg=0.620657
Tf=1.855542
Tg=0.620622
Tf=1.851491
Tg=0.620602
Tf=1.848343
Tg=0.620588
Tf=1.845805
Tg=0.620579
Tf=1.843703
Tg=0.620572
Tf=1.841925
Tg=0.620567
Tfbar=1.809048
Tgbar=0.620537
errof=
  列 1 至 2
   0.103950   0.073511
  列 3 至 4
   0.060023   0.051982
  列 5 至 6
   0.046494   0.042443
  列 7 至 8
   0.039295   0.036757
  列 9 至 10
   0.034655   0.032876
errog=
   1.0e-3 *
  列 1 至 2
   0.959757   0.339227
  列 3 至 4
   0.184628   0.119910
  列 5 至 6
   0.085796   0.065264
  列 7 至 8
   0.051790   0.042388
  列 9 至 10
   0.03552   0.030329

3. 变量代换后的复化梯形公式

编辑器:

syms x
f(x)=2*cos(x^2);
g(x)=2*sin(x^2);
errof=[];
errog=[];
for n=20:20:200
    h=1/n;
    x=linspace(1/n,1,n);
    Tf=eval((h*(0+sum(f(x))-f(x(1,n))/2)))
    Tg=eval((h*(0+sum(g(x))-g(x(1,n))/2)))
    errof(n/20)=abs(Tf-Tfbar);
    errog(n/20)=abs(Tg-Tgbar);
end
Tfbar=eval(int(f(x),0,1))
Tgbar=eval(int(g(x),0,1))
errof,errog

Tf=1.758347
Tg=0.620987
Tf=1.783873
Tg=0.620649
Tf=1.792304
Tg=0.620587
Tf=1.796505
Tg=0.620565
Tf=1.799020
Tg=0.620555
Tf=1.800696
Tg=0.620549
Tf=1.801891
Tg=0.620546
Tf=1.802785
Tg=0.620544
Tf=1.803484
Tg=0.620542
Tf=1.804041
Tg=0.620541
Tfbar=1.809048
Tgbar=0.620537
errof=
  列 1 至 2
   0.050701   0.025175
  列 3 至 4
   0.016745   0.012544
  列 5 至 6
   0.010028   0.008353
  列 7 至 8
   0.007157   0.006261
  列 9 至 10
   0.005564   0.005007
errog=
   1.0e-3 *
  列 1 至 2
   0.450502   0.112579
  列 3 至 4
   0.050031   0.028142
  列 5 至 6
   0.018010   0.012507
  列 7 至 8
   0.009189   0.007035
  列 9 至 10
   0.0055587   0.004503

4. Gauss-Legendre公式

编辑器：

syms x
f(x)=cos(((1+x)^2)/4);
g(x)=sin(((1+x)^2)/4);
errof=[];
errog=[];
x=[0,0,0,0,0;-0.5774,0.5774,0,0,0;0,-0.7746,0.7746,0,0;-0.8611,0.8611,-0.3400,0.3400,0;0,-0.9062,0.9062,-0.2309,0.2309,-0.5385,0.5385];
A=
[2,0,0,0,0;1,1,0,0,0;0.8889,0.5556,0.5556,0,0;0.3479,0.3479,0.6521,0.6521,0;0.5689,0.2369,0.2369,0.4786,0.4786];
for n=1:1:5
    S=0;
    R=0;
    for i=1:1:n
        S=S+A(n,i)*f(x(n,i));
        R=R+A(n,i)*g(x(n,i));
    end
    Tf=eval(S)
    Tg=eval(R)
    errof(n)=abs(Tf-Tfbar);
    errog(n)=abs(Tg-Tgbar);
end
Tfbar=eval(int(f(x),-1,1))
Tgbar=eval(int(g(x),-1,1))
errof,errog

Tf=1.937825
Tg=0.494808
Tf=1.811690
Tg=0.627333
Tf=1.808942
Tg=0.620590
Tf=1.809044
Tg=0.620540
Tf=1.808954
Tg=0.620509
Tfbar=1.809048
Tgbar=0.620537
errof=
  列 1 至 2
   0.128776   0.002642
  列 3 至 4
   0.000107   0.000001
  列 5
   0.000095
errog=
  列 1 至 2
   0.125729   0.006797
  列 3 至 4
   0.000053   0.000003
  列 5
   0.000027

5. Gauss-Jocabi公式

1. 正交多项式

syms x
w(x)=1/sqrt(x);
syms a [1,5]
syms P [1,5]
for n=1:1:5
    s=1;
    for i=1:1:n
        s=a(i)+s*x;
    end
    P(i)=s;
end
P

n=1
clear eq
syms eq [1,n]
outfor j=1:1:n
        eq(j)=int(w(x)*x^(j-1)*P(n),0,1);
    end
out=solve(eq);
out(1)

n=2
clear eq
syms eq [1,n]
for j=1:1:n
        eq(j)=int(w(x)*x^(j-1)*P(n),0,1);
    end
out=solve(eq);
out.a1,out.a2

n=3
clear eq
syms eq [1,n]
for j=1:1:n
        eq(j)=int(w(x)*x^(j-1)*P(n),0,1);
    end
out=solve(eq);
out.a1,out.a2,out.a3

n=4
clear eq
syms eq [1,n]
for j=1:1:n
        eq(j)=int(w(x)*x^(j-1)*P(n),0,1);
    end
out=solve(eq);
out.a1,out.a2,out.a3,out.a4

n=5
clear eq
syms eq [1,n]
for j=1:1:n
        eq(j)=int(w(x)*x^(j-1)*P(n),0,1);
    end
out=solve(eq);
out.a1,out.a2,out.a3,out.a4,out.a5

P = 
 [a1+x,a2+x*(a1+x),a3+x*(a2+x*(a1+x)),a4+x*(a3+x*(a2+x*(a1+x))),a5+x*(a4+x*(a3+x*(a2+x*(a1+x))))]
n =
    1
ans = 
   -1/3
n =
    2
ans =
   -6/7
ans =
   3/35
n =
   3
ans =
   -15/11
ans =
   5/11
ans =
   -5/231
n =
    4
ans = 
   -28/15
ans =
   14/13
ans =
   -28/143
ans =
   7/1287
n =
    5
ans = 
   -45/19
ans =
   630/323
ans =
   -210/323
ans =
   315/4199
ans =
   -63/46189

2. 零点

roots([1,-1/3])
roots([1,-6/7,3/35])
roots([1,-15/11,5/11,-5/231])
roots([1,-28/15,14/13,-28/143,7/1287])
roots([1,-45/19,630/323,-210/323,315/4199,-63/46189])

ans =  0.33333
ans =

   0.74156
   0.11559

ans =

   0.869499
   0.437198
   0.056939

ans =

   0.922157
   0.634677
   0.276184
   0.033648

ans =

   0.948494
   0.748335
   0.461597
   0.187832
   0.022164

3. 积分系数

X=[1/3,0,0,0,0;0.115587,0.741556,0,0,0;0.056939,0.437198,0.869500,0,0;0.033648,0.276184,0.634677,0.748335,0.948494];
syms x
w(x)=1/x^(1/2)
syms a [1,5]
syms P [1,5]
for n=1:1:5
    n
    s=1;
    for i=1:1:n
        s=s*(x-X(n,i));
    end
    w(x)=s;
    dw(x)=diff(w(x),1);
    for j=1:1:n
        j
        A=eval(int(w(x)*w(x)/((x-X(n,j))*dw(X(n,j))),0,1))
    end
end

n = 
    1
j = 
    1
A =
    2
n = 
    2
j = 
    1
A = 
    1.304290
j = 
    2
A = 
    0.695710
n = 
    3
j = 
    1
A = 
    0.935828
j = 
    2
A = 
    0.721523
j = 
    3
A = 
    0.342649
n = 
    4
j = 
    1
A =
    0.725368
j = 
    2
A = 
    0.627413
j = 
    3
A = 
    0.444762
j = 
    4
A = 
    0.202457
n = 
    5
j = 
    1
A = 
    0.591048
j = 
    2
A =
    0.538533
j = 
    3
A = 
    0.438173
j = 
    4
A = 
    0.298903 - 0.000000000000005i
j = 
    5
A = 
    0.133343 + 0.000000000000007i

4. 积分

syms x
f(x)=cos(x)
g(x)=sin(x)
F(x)=cos(x)/sqrt(x)
G(x)=sin(x)/sqrt(x)
X=[1/3,0,0,0,0;0.115587,0.741556,0,0,0;0.056939,0.437198,0.869499,0,0;0.033648,0.276184,0.634677,0.922157,0;0.022164,0.187832,0.461597,0.748335,0.948494];
A=
[2,0,0,0,0;1.304290,0.695710,0,0,0;0.935828,0.721523,0.342690,0,0;0.725368,0.627413,0.444762,0.202457,0;0.591048,0.538533,0.438127,0.298903,0.133343];
for n=1:1:5
    S=0;
    R=0;
    for i=1:1:n
        S=S+A(n,i)*f(X(n,i));
        R=R+A(n,i)*g(X(n,i));
    end
    Tf=eval(S)
    Tg=eval(R)
errof=[];
errog=[];
Tfbar=eval(int(F(x),0,1))
Tgbar=eval(int(G(x),0,1))
errof(n)=abs(Tf-Tfbar)
errog(n)=abs(Tg-Tgbar)
end 
errof,errog

X =
  列 1 至 2 
   0.333333 0 0.115587 0.741556 0.056939   0.437198 0.033648 0.276184 0.022164 0.187831     列 3 至 4 
   0 0 0 0 0.869500 0 0.634677 0.922157 0.461597 0.748335 
  列 5 
   0 0 0 0 0.948494 
A =
  列 1 至 2 
   2.000000 0 1.304290 0.695710 0.935828 0.721523 0.725368 0.627413 0.591048 0.538533   
  列 3 至 4 
   0 0 0 0 0.342649 0 0.444762 0.202457 0.438173 0.298903 
  列 5 
   0 0 0 0 0.133343 
Tf =
     1.889914
Tg =
     0.654389
Tf =
     1.808616
Tg =
     0.620331 
Tf =
     1.809049 
Tg =
     0.620537
Tf =
     1.809048 
Tg =
     0.620537 
Tf =
     1.809048 
Tg =
     0.620537 
errof = 
  列 1 至 2 
   0.080865 0.000432 
  列 3 至 4 
   0.00000091 0.00000000102
  列 5 
   0.000000000000659 
errog = 
  列 1 至 2 
   0.033852 0.000206 
  列 3 至 4 
   0.000000453 0.00000000052 
  列 5 
   0.000000000000295

4.结论

这五种方法插值积分的逼近精确度依次递增，所需迭代次数依次减少，收敛速度

逐渐加快，但计算复杂度也随之增加，且大多数Ic的逼近结果不如Is。

复化梯形公式，Newton-Cotes公式，变量代换后的复化梯形公式，Gauss-Legendre公式，Gauss-Jacobi公式插值积分的精确度比较

1.问题

2.数学理论和算法

1. n个等分区间的复化梯形公式n=100 : 100 : 1000

2. 带权的Newton-Cotes公式n=100 : 100 : ​1000

3. 变量代换后的n个等分区间的复化梯形公式n=20 : 20 : 200

4. Gauss-Legendre积分n=1 : 1 : 5

5. 由正交多项式推出相应Gauss型积分公式n=1 : 1 : 5

3.程序与结果

1. 复化梯形公式

编辑器：

2. 带权的Newton-Cotes公式

编辑器：

3. 变量代换后的复化梯形公式

编辑器:

4. Gauss-Legendre公式

编辑器：

5. Gauss-Jocabi公式

1. 正交多项式

2. 零点

3. 积分系数

4. 积分

4.结论

复化梯形公式，Newton-Cotes公式，变量代换后的复化梯形公式，Gauss-Legendre公式，Gauss-Jacobi公式插值积分的精确度比较

2. 带权的Newton-Cotes公式n=100 : 100 : 1000