梯形求积公式 和 复合梯形求积公式 Matlab 实现



梯形求积公式


  
仅使用区间两点x1,f(x1),x2,f(x2)

  组成的梯形面积S代替∫x2x1f(x)dx的近似方法 

∫x2x1f(x)dx≈S=x2−x12×(f(x2)−f(x1))



复合梯形求积公式


  
将求积区间[a,b]分为n个区间，每个区间步长为h(h=b−an)在每个区间求梯形积分 


  Si为第i个梯形的面积 
∫x2x1f(x)dx≈∑i=0nSi=h2∫(f(a)+∑k=1n−1f(xk)+f(b))



说明


  
上述公式是我用mathjax写的，如有错误请联系我修正  
敬请指正  

  概述省略了部分推导过程，请查阅详细推导资料




Matlab 实现代码



梯形求积公式


  
将该函数存为m文件




function res = Trapezium(f,a,b)
    format long;
    if b < a
        c = b;
        b = a;
        a = c;
    end
    res = (b-a)*(f(a)+f(b))/2;


  
调用下面语句测试函数




f = inline('sin(x)','x')
Trapezium(f,0,pi/2)



复合梯形求积公式


  
将该函数存为m文件




function res = ComTrapezium(f,n,a,b)
    format long;
    if b < a
        c = b;
        b = a;
        a = c;
    end
    h = (b-a)/n;
    d = f(a);
    for i = a+h:h:b-h 
        d = d + (2 * f(i));
    end
    d = d + f(b);
    res = (d * h / 2);


  
调用下面语句测试函数




f = inline('sin(x)','x')
ComTrapezium(f,4,0,pi/2)

一、简介

高斯求积公式是变步长数值积分的一种，基本形式是计算[-1,1]上的定积分。理论证明对于 n个节点的上述求积公式，最高有 2n - 1 次的代数精度，高斯公式就是使得上述公式具有 2n - 1次代数精度的积分公式。

详解



二、实现



# -*- coding: utf-8 -*-
"""
Created on Mon Dec 19 13:50:15 2016

    Gauss型求积公式

@author: Administrator
"""

from numpy import *
from scipy import integrate

import matplotlib.pyplot as plt

#(-1,1)
def f1(x):
    return x**2 / (1 - x**2)**(1/2)

#(0,pi/2]
def f2(x):
    return sin(x) / x

def f2_symbol(x):
    from sympy import sin
    return sin(x) / x

def gauss_legendre(f,a,b,n):
    from sympy import Symbol
    from sympy import diff
    from sympy import solve
    from sympy import Eq
    from math import factorial
    x = Symbol('x')
    l = diff((x**2 - 1)**(n+1),x,n+1)
    l = l  / (2**(n+1) * factorial(n+1))
    root = solve(Eq(l,0))
    A = []
    result = 0
    l_derivative = diff(l,x,1)
    for i in range(0,n+1):
        l_derivative_value = (l_derivative).evalf(subs={x:root[i]})
        A.append(2 / ((1 - root[i]**2) *  l_derivative_value**2))
        result += (b -a) / 2 * A[i] * f((a+b)/2 + (b-a)/2*root[i])
    return result

if __name__ == '__main__':
    n = 4 # 1 2 4
    print 'f1:'
    print 'Truth-value:',integrate.quad(f1,-1,1)[0]
    print 'Estimated-value:',float(gauss_legendre(f1,-1,1,n))
    print
    print 'f2:'
    print 'Truth-value',integrate.quad(f2,0,pi/2)[0]
    print 'Estimated-value:',float(gauss_legendre(f2_symbol,0,pi/2,n))
    zero = 1e-10
    x1 = linspace(-1,1,500)
    x1[0] += zero
    x1[-1] -= zero
    y1 = f1(x1)
    x2 = linspace(0,pi/2,500)
    x2[0] += zero
    y2 = f2(x2)
    fig = plt.figure(figsize=(8,6))
    ax1 = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax2.set_xlabel('x')
    ax2.set_ylabel('y')
    ax1.plot(x1,y1,color='b',linestyle='-',label='f1(x)')
    ax1.legend(loc='upper left')
    ax2.plot(x2,y2,color='r',linestyle='-',label='f2(x)')
    ax2.legend(loc='lower left')
    fig.show()
    fig.savefig('a.png')

为什么要引入数值积分
由于Gauss型求积公式属于数值积分的内容，学东西总要知道它的来龙去脉，下面我简单介绍一下为什么要引入数值积分
给定函数$f(x)\in C[a,b]$,考虑积分
$I(f)=\int_{a}^{b} f(x) dx$

的计算问题，从数学分析中知道，当已知$f(x)$的原函数为$F(x)$时，由牛顿-莱布尼兹公式，有

$I(f)=\int_{a}^{b}f(x)dx=F(b)-F(a)$

然而，在实际计算中，被积函数的$f(x)$的原函数经常无法用初等函数表示，但过于复杂。还有时，$f(x)$只在一些离散点上给出。在这样的情况下，就有必要借助数值方法来求$I(f)$的近似值。
Gauss型求积公式的定义
从正交多项式理论可知，在区间[a,b]上，对给定的权函数$\rho(x)$,存在正交多项式系${\omega_{k}(x)}_{k=0}^{\infty}$并且可以将其构造出来。进一步，已经证明了$\omega_{k}(x)$在[a,b]上恰有k个相异的根，取$\omega_{n+1}(x)$的n+1个根为求积结点，构造插值型求积公式，将得到具有2n+1次代数精度的求积公式，称如此构造出来的求积公式为Gauss型求积公式即：

如果求积结点$x_{0},x_{1},....,x_{n}$,使插值型求积公式$\int_{-1}^{1}f(x)dx\approx\sum_{k=0}^{n}A_kf(x_k)$的代数精度为2n+1，则称该求积公式为Gauss型求积公式，称这些求积结点为Gauss点。

特别地，Gauss求积公式当[a,b]=[-1,1]时，取$\rho(x)$=1,Gauss型求积公式称为Gauss求积公式，此时的求积结点

$x_k(k=1,2,3....)$为n次Legendre多项式的根。

且$A_k=\int_{-1}^{1} \frac{\omega(x)}{(x-x_k)\omega'(x_k)}dx$

而Legendre多项式为：



打代码太南了，我手写偷下懒。

Gauss公式的应用
Gauss公式主要应用为二点Gauss公式和三点Gauss公式

二点Gauss公式：当n=2时，有 $x_1=-\sqrt\frac{1}{3},x_2=\sqrt\frac{1}{3},A_1=A_2=1$,求积公式为$G_2(f)=f(-\sqrt\frac{1}{3})+f(\sqrt\frac{1}{3})$

三点Gauss求积公式：当n=3时，有$x_1=-\sqrt\frac{3}{5},x_2=0，x_3=\sqrt\frac{3}{5},A_1=A_3=\frac{5}{9},A_2=\frac{8}{9}$,求积公式为$G_3(f)=\frac{5}{9}f(-\sqrt\frac{3}{3})+\frac{8}{9}f(0)+\frac{5}{9}f(\sqrt\frac{3}{5})$

为了计算[a,b]上的定积分$\int_{a}^{b} f(x) dx$可以通过自变量的变换$x=\frac{1}{2}(a+b+t*(b-a))$将其转换为[-1,1]上的积分，然后使用Gauss积分公式计算其近似值。

例：

用两点Gauss公式计算$\int_{0}^{1}\frac{sinx}{x}dx$

解：变换令x=0.5(t+1)

$\int_{0}^{1}\frac{sinx}{x}dx=\int_{-1}^{1}\frac{sin0.5(t+1)}{t+1}dx$

取$t_0=\frac{-1}{\sqrt 3},t_1=\frac{1}{\sqrt 3}$

$\int_{0}^{1}\frac{sinx}{x}dx=0.5[\frac{sin0.5(t_0+1)}{0.5(t_0+1}+\frac{sin0.5(t_1+1)}{0.5(t_1+1}]$
Matlab程序
function [w,p] = Gauss_point_1D(n,a,b)
% Gauss quarature point on [-1,1]
if n == 2
    w = [1,1];
    p = [-1/sqrt(3),1/sqrt(3)];
elseif n == 4
     w = [0.3478548451,0.3478548451,0.6521451549,0.6521451549];
     p = [0.8611363116,-0.8611363116,0.3399810436,-0.3399810436];
elseif n == 8
     w = [0.1012285363,0.1012285363,0.2223810345,0.2223810345,0.3137066459,0.3137066459,0.3626837834,0.3626837834];
     p = [0.9602898565,-0.9602898565,0.7966664774,-0.7966664774,0.5255324099,-0.5255324099,0.1834346425,-0.1834346425];
end

% Gauss quarature point on [a,b]
w = 0.5*(b-a)*w;
p = 0.5*(b-a)*p+0.5*(b+a);



function q = Gauss_quadrature(fun,a,b,n)
% Gauss quadrature
[w,p] = Gauss_point_1D(n,a,b);
q = sum(w.*fun(p));

Gauss型求积公式
若机械求积公式具有阶代数精度，则称为Gauss型求积公式，而在上关于权函数的次正交多项式的零点就是Gauss型求积公式的Gauss点。
在Gauss型求积公式中，若权函数，区间为，则公式为
                
特别的称为Gauss-Legendre公式。下表列出Gauss-Legendre公式的结点和系数。
 
 
 
 
 
 

0


0.0000000


2.0000000


3


±0.8611361
±0.3399810


0.3478548
0.6521452


1


±0.5773503


1.000000


4


±0.9061798
±0.5384693
0.000000


0.2369269
0.4786287
0.5688889


2


±0.7745967
0.000000


0.5555556
0.8888889


5


±0.9324695
±0.6612094
±0.2386192


0.1713245
0.3607618
0.4679139

当积分区间是一般的区间时，只要做变换
                    
可将公式转换为
               
对等式右端积分使用Gauss-Legendre公式即可。
     Matlab代码如下：
function I = IntGaussLegen(f,a,b,n,AK,XK)
%
syms t;
t= findsym(sym(f));
if(n<5 && nargin == 4)
    AK = 0;
    XK = 0;
else                       
    XK1=((b-a)/2)*XK+((a+b)/2);
    I=((b-a)/2)*sum(AK.*subs(sym(f),findsym(f),XK1));
end
 
ta = (b-a)/2;
tb = (a+b)/2;
switch n
    case 0,
        I=2*ta*subs(sym(f),t,tb);
       
    case 1,
        I=ta*(subs(sym(f),t,ta*0.5773503+tb)+...
            subs(sym(f),t,-ta*0.5773503+tb));
       
    case 2,
        I=ta*(0.55555556*subs(sym(f),t,ta*0.7745967+tb)+...
            0.55555556*subs(sym(f),t,-ta*0.7745967+tb)+...
            0.88888889*subs(sym(f),t,tb));
          
    case 3,
        I=ta*(0.3478548*subs(sym(f),t,ta*0.8611363+tb)+...
            0.3478548*subs(sym(f),t,-ta*0.8611363+tb)+...
            0.6521452*subs(sym(f),t,ta*0.3398810+tb) +...
            0.6521452*subs(sym(f),t,-ta*0.3398810+tb));
         
    case 4,
        I=ta*(0.2369269*subs(sym(f),t,ta*0.9061793+tb)+...
            0.2369269*subs(sym(f),t,-ta*0.9061793+tb)+...
            0.4786287*subs(sym(f),t,ta*0.5384693+tb) +...
            0.4786287*subs(sym(f),t,-ta*0.5384693+tb)+...
            0.5688889*subs(sym(f),t,tb));
case 5,
        I=ta*(0.1713245*subs(sym(f),t,ta*0.9324695+tb)+...
            0.1713245*subs(sym(f),t,-ta*0.9324695+tb)+...
            0.3607616*subs(sym(f),t,ta*0.6612094+tb)+...
            0.3607616*subs(sym(f),t,-ta*0.6612094+tb)+...
            0.4679139*subs(sym(f),t,ta*0.2386292+tb)+...
            0.4679139*subs(sym(f),t,-*0.2386292+tb));
end
     若求积区间为，权函数为，则所建立的Gauss型求积公式为
称为Gauss-Chebyshev公式。其中结点为次Chebyshev多项式的零点，公式可记为
                    
 
    若求积区间为，权函数为，则所建立的Gauss型求积公式为
称为Gauss-Lagurre公式。其中结点为次Lagurre多项式的零点，下表列出Gauss-Lagurre公式的结点和系数。
 
 
 
 
 
 

0


1.0000000


1.0000000


 
 


0.2635603
1.4134031


0.52175556
0.3986668


1


0.5857864
3.4142135


0.8535533
0.1464466


4


3.5964258
7.0858100
12.6408008


0.0759425
0.00361176
0.00002337


2


0.4157746
2.2942804
6.2899451


0.7110930
0.2785177
0.0103893


5


0.2228466
1.1889321
2.9927363


0.4589647
0.4170008
0.1133734


3


0.3225477
1.7457611
4.5366203
9.3950710


0.6031541
0.3574187
0.0388879
0.0005393


 


5.7751436
9.8374674
15.9828740


0.1039920
0.000261017
0.0000008985

 
 
转载于:https://www.cnblogs.com/wander-clouds/p/9949154.html