1.问题

分别计算积分
$$I_c={\int }_0^1\frac{\cos x}{\sqrt{x}}dx=1.809048475800...$$

$$I_s={\int }_0^1\frac{\sin x}{\sqrt{x}}dx=0.620536603446$$

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

2.数学理论和算法

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

$$T_n=\sum _{k=0}^{n-1}\frac{(x_{k+1}-x_k)}{2}(f(x_k)+f(x_{k+1}))$$

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

推导过程：
$$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+ih,i=0,1,\cdot \cdot \cdot ，n,h=\frac{b-a}{n}$$

$$P_n(x)=\sum _{i=0}^n\frac{\omega (x)}{(x-x_i){\omega }^{\prime }(x)}f(x_i),\omega (x)=(x-x_0)(x-x_1)\cdot \cdot \cdot (x-x_n)$$

$${\int }_a^{b}f(x)dx\approx {\int }_a^b{P}_n(x)dx={\int }_a^b(\sum _{i=0}^n\frac{\omega (x)}{(x-x_i){\omega }^{\prime }(x_i)}f(x_i))dx=\sum _{i=0}^n({\int }_a^b\frac{\omega (x)}{(x-x_i){\omega }^{\prime }(x_i)}dx)f(x_i)=\sum _{i=0}^n{A}_{i}f(x_i)$$

$$A_i={\int }_a^b\frac{\omega (x)}{(x-x_i){\omega }^{\prime }(x_i)}dx$$

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

$$x=t^2,I_c^{\prime }={\int }_0^{1}2\cos (t^2)dt,I_s^{\prime }={\int }_0^{1}2\sin (t^2)dt$$

复化梯形公式如1.

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

$$t=\frac{0+1}{2}+\frac{1-0}{2}T=\frac{1+T}{2}$$

$$I_c={\int }_{-1}^1\cos (\frac{1+T}{2}{)}^{2}dT$$

$$I_s={\int }_{-1}^1\sin (\frac{1+T}{2}{)}^{2}dT$$

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:1000h=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 至 40.084315   0.073018列 5 至 60.065309   0.059619列 7 至 80.055196   0.051631列 9 至 100.048679   0.046181
errog=1.0e-03 *列 1 至 20.206890   0.073250列 3 至 40.039897   0.025923列 5 至 60.018554   0.014117列 7 至 80.011204   0.009172列 9 至 100.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:1000h=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 至 20.103950   0.073511列 3 至 40.060023   0.051982列 5 至 60.046494   0.042443列 7 至 80.039295   0.036757列 9 至 100.034655   0.032876
errog=1.0e-3 *列 1 至 20.959757   0.339227列 3 至 40.184628   0.119910列 5 至 60.085796   0.065264列 7 至 80.051790   0.042388列 9 至 100.03552   0.030329

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

编辑器:

syms x
f(x)=2*cos(x^2);
g(x)=2*sin(x^2);
errof=[];
errog=[];
for n=20:20:200h=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 至 20.050701   0.025175列 3 至 40.016745   0.012544列 5 至 60.010028   0.008353列 7 至 80.007157   0.006261列 9 至 100.005564   0.005007
errog=1.0e-3 *列 1 至 20.450502   0.112579列 3 至 40.050031   0.028142列 5 至 60.018010   0.012507列 7 至 80.009189   0.007035列 9 至 100.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:5S=0;R=0;for i=1:1:nS=S+A(n,i)*f(x(n,i));R=R+A(n,i)*g(x(n,i));endTf=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 至 20.128776   0.002642列 3 至 40.000107   0.000001列 50.000095
errog=列 1 至 20.125729   0.006797列 3 至 40.000053   0.000003列 50.000027

5. Gauss-Jocabi公式

1. 正交多项式

$$

$$

syms x
w(x)=1/sqrt(x);
syms a [1,5]
syms P [1,5]
for n=1:1:5s=1;for i=1:1:ns=a(i)+s*x;endP(i)=s;
end
Pn=1
clear eq
syms eq [1,n]
outfor j=1:1:neq(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:neq(j)=int(w(x)*x^(j-1)*P(n),0,1);end
out=solve(eq);
out.a1,out.a2n=3
clear eq
syms eq [1,n]
for j=1:1:neq(j)=int(w(x)*x^(j-1)*P(n),0,1);end
out=solve(eq);
out.a1,out.a2,out.a3n=4
clear eq
syms eq [1,n]
for j=1:1:neq(j)=int(w(x)*x^(j-1)*P(n),0,1);end
out=solve(eq);
out.a1,out.a2,out.a3,out.a4n=5
clear eq
syms eq [1,n]
for j=1:1:neq(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.741560.11559ans =0.8694990.4371980.056939ans =0.9221570.6346770.2761840.033648ans =0.9484940.7483350.4615970.1878320.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:5ns=1;for i=1:1:ns=s*(x-X(n,i));endw(x)=s;dw(x)=diff(w(x),1);for j=1:1:njA=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:5S=0;R=0;for i=1:1:nS=S+A(n,i)*f(X(n,i));R=R+A(n,i)*g(X(n,i));endTf=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。