【数值分析】高斯型求积公式,任意区间三点gauss求积公式,matlab实现

你哥同学 2024-07-02 11:35:03 阅读 53

Gauss型求积公式

Gauss型求积公式定义

a

b

ρ

(

x

)

f

(

x

)

d

x

i

=

1

n

A

i

f

(

x

i

)

\int_{ a }^{b} \rho(x)f(x) \mathrm dx \approx \sum_{i=1}^{ n}A_if(x_i)

∫ab​ρ(x)f(x)dx≈i=1∑n​Ai​f(xi​)

如果求积公式具有

2

n

1

{2n-1}

2n−1 次代数精度,则称对应的节点

x

1

,

x

2

,

,

x

n

{x_1,x_2,\cdots ,x_n}

x1​,x2​,⋯,xn​ 为Gauss点,此时求积公式称为Gauss型求积公式。

为了讨论方便,本节取

n

{n}

n 个节点,并记节点为

x

1

,

x

2

,

x

n

{x_1,x_2,\cdots x_n}

x1​,x2​,⋯xn​ ,从

1

{1}

1 开始取!同时,所讨论的积分均为带有权函数

ρ

(

x

)

{\rho(x)}

ρ(x) 的积分。

插值型求积公式

a

b

ρ

(

x

)

f

(

x

)

d

x

i

=

1

n

A

i

f

(

x

i

)

\int_{ a }^{b} \rho(x)f(x) \mathrm dx \approx \sum_{i=1}^{ n}A_if(x_i)

∫ab​ρ(x)f(x)dx≈i=1∑n​Ai​f(xi​)

的代数精度最高不超过

2

n

1

{2n-1}

2n−1 ,且达到

2

n

1

{2n-1}

2n−1 时,所有求积系数为正。所以高斯公式是给定节点数下代数精度最高的求积公式。

构造Gauss型求积公式的步骤

对给定区间

[

a

,

b

]

{[a,b]}

[a,b] 及权函数

ρ

(

x

)

{\rho(x)}

ρ(x) ,由Schmidt正交化过程构造正交多项式

P

0

(

x

)

,

P

1

(

x

)

,

,

P

n

(

x

)

{P_0(x),P_1(x), \cdots ,P_n(x)}

P0​(x),P1​(x),⋯,Pn​(x)求出

P

n

(

x

)

{P_n(x)}

Pn​(x) 的

n

{n}

n 个零点

x

1

,

x

2

,

,

x

n

{ x_1,x_2, \cdots, x_n}

x1​,x2​,⋯,xn​ 即为Gauss点计算求积系数

A

i

=

a

b

ρ

(

x

)

l

i

(

x

)

d

x

{A_i= \int_{ a }^{b} \rho(x)l_i(x) \mathrm dx }

Ai​=∫ab​ρ(x)li​(x)dx,

i

=

1

,

2

,

,

n

{ i=1,2,\cdots,n }

i=1,2,⋯,n,

l

i

{l_i}

li​ 为拉格朗日基函数

[!example]-

求计算积分

1

1

x

2

f

(

x

)

d

x

{ \int_{ -1 }^{1} x^2f(x) \mathrm dx}

∫−11​x2f(x)dx 的两点Gauss公式。

解:

ρ

(

x

)

=

x

2

,

n

=

2

\rho(x)=x^2, n=2

ρ(x)=x2,n=2

首先按Schmidt正交化求出正交多项式:

P

0

(

x

)

=

1

P

1

(

x

)

=

x

(

x

,

P

0

(

x

)

)

(

P

0

(

x

)

,

P

0

(

x

)

)

P

0

(

x

)

=

x

1

1

x

3

d

x

1

1

x

2

d

x

=

x

P

2

(

x

)

=

x

2

(

x

2

,

P

0

(

x

)

)

(

P

0

(

x

)

,

P

0

(

x

)

)

P

0

(

x

)

(

x

2

,

P

1

(

x

)

)

(

P

1

(

x

)

,

P

1

(

x

)

)

P

1

(

x

)

=

x

2

1

1

x

4

d

x

1

1

x

2

d

x

1

1

x

5

d

x

1

1

x

4

d

x

x

=

x

2

3

5

\begin{align*} P_0(x)=&1 \\ \\ P_1(x)=&x- \frac{(x,P_0(x))}{(P_0(x),P_0(x))}P_0(x)=x- \frac{\int_{ -1 }^{1} x^3 \mathrm dx}{\int_{ -1 }^{1} x^2 \mathrm dx}=x \\ \\ P_2(x)=&x^2- \frac{(x^2,P_0(x))}{(P_0(x),P_0(x))}P_0(x)-\frac{(x^2,P_1(x))}{(P_1(x),P_1(x))}P_1(x) \\ \\ =&x^2- \frac{\int_{ -1 }^{1} x^4 \mathrm dx}{\int_{ -1 }^{1} x^2 \mathrm dx}- \frac{\int_{ -1 }^{1} x^5 \mathrm dx}{\int_{ -1 }^{1} x^4 \mathrm dx}x=x^2- \frac{3}{5} \end{align*}

P0​(x)=P1​(x)=P2​(x)==​1x−(P0​(x),P0​(x))(x,P0​(x))​P0​(x)=x−∫−11​x2dx∫−11​x3dx​=xx2−(P0​(x),P0​(x))(x2,P0​(x))​P0​(x)−(P1​(x),P1​(x))(x2,P1​(x))​P1​(x)x2−∫−11​x2dx∫−11​x4dx​−∫−11​x4dx∫−11​x5dx​x=x2−53​​

再令

P

2

(

x

)

=

0

{P_2(x)=0}

P2​(x)=0 求出Gauss点:

x

1

=

3

5

,

x

2

=

3

5

x_1=-\sqrt{\frac{3}{5}} \,\,,\,\, x_2=\sqrt{\frac{3}{5}}

x1​=−53​

​,x2​=53​

最后计算求积系数:

A

1

=

1

1

x

2

l

1

(

x

)

d

x

=

1

1

x

2

x

x

2

x

1

x

2

d

x

=

1

3

A

2

=

1

1

x

2

l

2

(

X

)

d

x

=

1

1

x

2

x

x

1

x

1

x

2

d

x

=

1

3

\begin{align*} A_1=& \int_{ -1 }^{1} x^2l_1(x) \mathrm dx= \int_{ -1 }^{1} x^2 \frac{x-x_2}{x_1-x_2} \mathrm dx= \frac{1}{3}\\ \\ A_2=& \int_{ -1 }^{1} x^2l_2(X) \mathrm dx= \int_{ -1 }^{1} x^2 \frac{x-x_1}{x_1-x_2} \mathrm dx= \frac{1}{3} \end{align*}

A1​=A2​=​∫−11​x2l1​(x)dx=∫−11​x2x1​−x2​x−x2​​dx=31​∫−11​x2l2​(X)dx=∫−11​x2x1​−x2​x−x1​​dx=31​​

所以两点Gauss公式为

1

1

x

2

f

(

x

)

d

x

1

3

[

f

(

3

5

)

+

f

(

3

5

)

]

\int_{ -1 }^{1} x^2f(x) \mathrm dx \approx \frac{1}{3}[f(- \sqrt{\frac{3}{5}})+f(\sqrt{\frac{3}{5}})]

∫−11​x2f(x)dx≈31​[f(−53​

​)+f(53​

​)]

4.1 Gauss-legendre勒让德求积公式

区间

[

1

,

1

]

{[-1,1]}

[−1,1] ,权函数

ρ

(

x

)

=

1

{\rho(x)=1}

ρ(x)=1 。

Gauss求积公式中在区间

[

1

,

1

]

{[-1,1]}

[−1,1] 上权函数为

1

{1}

1 的Gauss点其实是确定的,所以可以通过把待求区间伸缩变换到

[

1

,

1

]

{[-1,1]}

[−1,1] 间,再通过Gauss求积公式求解。

一点Gauss求积公式

1

1

f

(

x

)

d

x

2

f

(

0

)

\int_{ -1 }^{1} f(x) \mathrm dx \approx 2f(0)

∫−11​f(x)dx≈2f(0)

两点Gauss求积公式

1

1

f

(

x

)

d

x

f

(

1

3

)

+

f

(

1

3

)

\int_{ -1 }^{1} f(x) \mathrm dx \approx f(- \frac{1}{\sqrt{3}})+f(\frac{1}{\sqrt{3}})

∫−11​f(x)dx≈f(−3

​1​)+f(3

​1​)

三点Gauss求积公式

1

1

f

(

x

)

d

x

5

9

f

(

3

5

)

+

8

9

f

(

0

)

+

5

9

f

(

3

5

)

\int_{ -1 }^{1} f(x) \mathrm dx \approx \frac{5}{9} f(- \sqrt{\frac{3}{5}})+ \frac{8}{9}f(0) +\frac{5}{9}f(\sqrt{\frac{3}{5}})

∫−11​f(x)dx≈95​f(−53​

​)+98​f(0)+95​f(53​

​)

任意区间的三点高斯求积公式:

a

b

f

(

x

)

d

x

(

b

a

)

2

(

5

9

f

(

(

a

+

b

)

2

3

5

(

b

a

)

2

)

+

8

9

f

(

(

a

+

b

)

2

)

+

5

9

f

(

(

a

+

b

)

2

+

3

5

(

b

a

)

2

)

)

\int_{ a }^{b} f(x) \mathrm dx \approx \frac{(b-a)}{2} \bigg( \frac{5}{9} f(\frac{(a+b)}{2} - \sqrt{\frac{3}{5}}\frac{(b-a)}{2})+ \frac{8}{9}f(\frac{(a+b)}{2}) +\frac{5}{9}f(\frac{(a+b)}{2} + \sqrt{\frac{3}{5}}\frac{(b-a)}{2}) \bigg)

∫ab​f(x)dx≈2(b−a)​(95​f(2(a+b)​−53​

​2(b−a)​)+98​f(2(a+b)​)+95​f(2(a+b)​+53​

​2(b−a)​))

matlab实现

%% 三点高斯勒让德求积公式

% 输入函数,积分上界,积分下界

function I = gaussL3P(f,a,b)

x = (a+b)/2+(b-a)/2*[-sqrt(0.6) 0 sqrt(0.6)];

I = (b-a)/2*f(x)*[5 8 5]'/9;

end

四点Gauss求积公式

高斯点

:

±

0.8611363

,

高斯系数

:

0.3478548

高斯点:\pm0.8611363 \,\,,\,\, 高斯系数:0.3478548

高斯点:±0.8611363,高斯系数:0.3478548

高斯点

:

±

0.3399810

,

高斯系数

:

0.6521452

高斯点:\pm0.3399810 \,\,,\,\, 高斯系数:0.6521452

高斯点:±0.3399810,高斯系数:0.6521452

[!example]-

已知三点Gauss公式

1

1

f

(

x

)

d

x

5

9

f

(

0.6

)

+

8

9

f

(

0

)

+

5

9

f

(

0.6

)

\int_{ -1 }^{1} f(x) \mathrm dx \approx \frac{5}{9}f(\sqrt{0.6})+ \frac{8}{9}f(0)+ \frac{5}{9}f(- \sqrt{0.6})

∫−11​f(x)dx≈95​f(0.6

​)+98​f(0)+95​f(−0.6

​)

试用如上公式计算

0.5

1

x

d

x

{ \int_{ 0.5 }^{1} \sqrt{x} \mathrm dx}

∫0.51​x

​dx 的值。

解:换限,设

t

=

a

x

+

b

t=ax+b

t=ax+b

{

0.5

a

+

b

=

1

a

+

b

=

1

{

a

=

4

b

=

3

\begin{cases} 0.5a+b=-1 \\ \\ a+b=1 \end{cases} \Rightarrow \begin{cases} a=4 \\ \\ b=-3 \end{cases}

⎧​0.5a+b=−1a+b=1​⇒⎩

⎧​a=4b=−3​

0.5

1

f

(

x

)

d

x

=

1

4

1

1

f

(

t

+

3

4

)

d

t

1

4

[

5

9

f

(

0.6

+

3

4

)

+

8

9

f

(

0

+

3

4

)

+

5

9

f

(

0.6

+

3

4

)

]

=

1

4

[

5

9

0.6

+

3

4

+

8

9

0

+

3

4

+

5

9

0.6

+

3

4

]

=

0.4310

\begin{align*} \int_{ 0.5 }^{1} f(x) \mathrm dx=& \frac{1}{4} \int_{ -1 }^{1} f(\frac{t+3}{4}) \mathrm dt \\ \\ \approx& \frac{1}{4}[\frac{5}{9} f(\frac{\sqrt{0.6}+3}{4})+ \frac{8}{9}f(\frac{0+3}{4}) + \frac{5}{9}f(\frac{- \sqrt{0.6}+3}{4}) ] \\ \\ =&\frac{1}{4}[\frac{5}{9} \sqrt{\frac{\sqrt{0.6}+3}{4}}+ \frac{8}{9}\sqrt{\frac{0+3}{4}} + \frac{5}{9}\sqrt{\frac{- \sqrt{0.6}+3}{4}} ] \\ \\ =&0.4310 \end{align*}

∫0.51​f(x)dx=≈==​41​∫−11​f(4t+3​)dt41​[95​f(40.6

​+3​)+98​f(40+3​)+95​f(4−0.6

​+3​)]41​[95​40.6

​+3​

​+98​40+3​

​+95​4−0.6

​+3​

​]0.4310​

4.2 变形的Gauss型求积公式

Gauss-Laguerre拉盖尔求积公式

区间

[

0

,

+

)

{[0,+\infty)}

[0,+∞) ,权函数

ρ

(

x

)

=

e

x

{\rho(x)=e^{-x}}

ρ(x)=e−x

Gauss-Hermite求积公式

两点Gauss-Hermite求积公式

+

e

x

2

f

(

x

)

d

x

π

2

f

(

2

2

)

+

π

2

f

(

2

2

)

\int_{ -\infty }^{+\infty} e^{-x^2}f(x) \mathrm dx \approx \frac{\sqrt{\pi}}{2} f(-\frac{\sqrt{2}}{2})+\frac{\sqrt{\pi}}{2} f(\frac{\sqrt{2}}{2})

∫−∞+∞​e−x2f(x)dx≈2π

​​f(−22

​​)+2π

​​f(22

​​)

三点Gauss-Hermite求积公式

+

e

x

2

f

(

x

)

d

x

π

6

f

(

6

2

)

+

2

π

3

f

(

0

)

+

π

6

f

(

6

2

)

\int_{ -\infty }^{+\infty} e^{-x^2}f(x) \mathrm dx \approx \frac{\sqrt{\pi}}{6} f(-\frac{\sqrt{6}}{2})+ \frac{2 \sqrt{\pi}}{3} f(0)+\frac{\sqrt{\pi}}{6} f(\frac{\sqrt{6}}{2})

∫−∞+∞​e−x2f(x)dx≈6π

​​f(−26

​​)+32π

​​f(0)+6π

​​f(26

​​)

区间

(

,

+

)

{(-\infty,+\infty)}

(−∞,+∞) ,权函数

ρ

(

x

)

=

e

x

2

{\rho(x)=e^{-x^2}}

ρ(x)=e−x2



声明

本文内容仅代表作者观点,或转载于其他网站,本站不以此文作为商业用途
如有涉及侵权,请联系本站进行删除
转载本站原创文章,请注明来源及作者。