精华内容
下载资源
问答
  • 梯形求积公式 和 复合梯形求积公式 Matlab 实现梯形求积公式 仅使用区间两点x1,f(x1),x2,f(x2)x_1,f(x_1),x_2,f(x_2) 组成的梯形面积S代替∫x2x1f(x)dx \int_{x_1}^{x_2} f(x)dx 的近似方法 ∫x2x1f(x)dx≈S=...

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

    梯形求积公式

    仅使用区间两点x1,f(x1),x2,f(x2)
    组成的梯形面积S代替x2x1f(x)dx的近似方法
    求积

    x2x1f(x)dxS=x2x12×(f(x2)f(x1))

    复合梯形求积公式

    将求积区间[a,b]分为n个区间,每个区间步长为h(h=ban)在每个区间求梯形积分
    复合梯形求积
    Si为第i个梯形的面积

    x2x1f(x)dxi=0nSi=h2(f(a)+k=1n1f(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)
    展开全文
  • matlab 复合梯形求积公式,.M文件,可直接运行出结果。
  • 【C语言基础】利用复合梯形求积公式计算定积分 一、复合梯形求积公式 这是数值分析中一种求解定积分的近似方法。适用于被积函数的原函数不能用初等函数表示的情况。 基本思路 将被积函数 f(x)与x轴围成的区域分成n...

    【C语言基础】利用复合梯形求积公式计算定积分

    一、复合梯形求积公式

    这是数值分析中一种求解定积分的近似方法。适用于被积函数的原函数不能用初等函数表示的情况。

    基本思路

    将被积函数 f(x)与x轴围成的区域分成n个梯形,把n个梯形面积求和得到积分的近似值。若精度不满足需要,则可以将每个区间再等分一次,得到2n+1个等分区间,然后再求和,直到精度满足需要。

    公式

    abf(x)dxh2[f(a)+f(b)+2i=1n1f(xi)] \int\limits_a^b f(x)dx \approx \frac{h}{2}[f(a)+f(b)+2\sum_{i=1}^{n-1} f(x_i)]
    其中:

    n 为等分区间数
    h=(b-a)/2 (积分步长)

    C语言实现算法

    以下面积分为例:

    011+exsin4xdx \int\limits_0^11+e^-xsin4xdx

    首先写出计算被积函数值的函数
    double f(double x)
    {
        return 1 + exp(-x)*sin(4x);
    }
    
    然后写出求积分近似解的表达式
    double caculate(double a,double b)
    {
        double y1,y2,h;
        int n=1,i;
        h = a-b;
        y1 = h*(f(a)+f(b))/2;
        do
        {
           y2 = y1;
           n = n*2;
           h = (a-b)/n;
           y1 = f(a)+f(b);
           for(i=1;i<n:i++)
                 y1 = y1 + f(a + i*h);
            y1 = h*y1/2;
        }while(fabs(y2-y1)>1.0e-5));
        return y2;
    }
    
    其中frab()是用来计算绝对值的函数
    展开全文
  • 复合梯形公式和复合辛普森求积公式计算。      ​ 复合梯度代码 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function T_n=fht(a,b,n) h=(b-a)/n; for k=0:n x(k+1)=a+k*h; if...
    • 实验内容

    1. 计算圆周率π,要求误差小于10-8。
    2. 用复合梯形公式和复合辛普森求积公式计算。  

     

     

    ​
    
    
    复合梯度代码
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function T_n=fht(a,b,n) 
    h=(b-a)/n; 
    for k=0:n     
        x(k+1)=a+k*h;     
        if x(k+1)==0          
            x(k+1)=10^(-10);     
        end
    end
    T_1=h/2*(fx1(x(1))+fx1(x(n+1))); 
    for i=2:n      
        F(i)=h*fx1(x(i)); 
    end
    T_2=sum(F); 
    T_n=T_1+T_2;
    
      复合辛普森求积代码
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function S_n=S_P_S(a,b,n) 
    h=(b-a)/n;
    for k=0:n     
        x(k+1)=a+k*h;      
        x_k(k+1)=x(k+1)+1/2*h;     
        if (x(k+1)==0)|(x_k(k+1)==0)         
            x(k+1)=10^(-10);         
            x_k(k+1)=10^(-10); 
        end
    end
    S_1=h/6*(fx1(x(1))+fx1(x(n+1)));
    for i=2:n      
         F_1(i)=h/3*fx1(x(i)); 
    end
    for j=1:n      
         F_2(j)=2*h/3*fx1(x_k(j)); 
    end
    S_2=sum(F_1)+sum(F_2); 
    S_n=S_1+S_2;
                                         龙贝格和变步长公式代码
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    function asd
    format long;
       a=0;
       b=1;
       m=2;
       t(1)=0.5*(b-a)*(f(a)+f(b));
       t(2)=0.5*t(1)+0.5*(b-a)*f((a+b)/2);
       s(1)=(4*t(2)-t(1))/3;
       j=2;
       while abs (t(j)-t(j-1))>(0.5e-8/4),
           h=(b-a)/m;
           k=0:(m-1);
           j=j+1;
           t(j)=0.5*t(j-1)+0.5*h*sum(f(a+(k+1/2)*h));
           s(j-1)=(4*t(j)-t(j-1))/3;
           c(j-2)=(16*s(j-1)-s(j-2))/15;
           if j>3,
               r(j-3)=(64*c(j-2)-c(j-3))/63;
           end
           m=m*2;
       end
       fsimpson=4*s(end)
       rbg=4*r(end)
       bbb=4*c(end)
       j
                                      方程代码
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    function y=fx1(x) 
    y=4/(1+x^2);
    
    [点击并拖拽以移动]
    ​

     

    >> T_20=fht(0,1,1048576)
    
    T_20 =
    
        3.1416
    
    >> vpa(T_20,20)
     
    ans =
     
    3.1415926535896492311
     
    >> 
    
    
    
    
    S_1=S_P_S(0,1,2)
    
    S_1 =
    
        3.2200
    
    >>  vpa(S_1,7) 
     
    ans =
     
    3.22
     
    >>  S_2=S_P_S(0,1,4)
    
    S_2 =
    
        3.1518
    
    >> vpa(S_2,7)
     
    ans =
     
    3.151849
     
    >> S_3=S_P_S(0,1,8)
    
    S_3 =
    
        3.1429
    
    >> vpa(S_3,7)
     
    ans =
     
    3.14289
     
    >> S_4=S_P_S(0,1,16)
    
    S_4 =
    
        3.1418
    
    
    bbb =
    
       3.1415955553589793
    
    
    bbb =
    
       3.14159244666653589793
    
    bbb =
    
       3.141592653589793

     

    展开全文
  • 1.求积公式余项 1.1 定义 R[f]=∫ab ⁣ ⁣ ⁣f(x)dx−∑k=0nAkf(xk)=Kf(m+1)(η),(1)R[f]=\int_a^b\!\!\!f(x)\mathrm{d}x-\sum_{k=0}^nA_kf(x_k)=Kf^{(m+1)}(\eta ),(1)R[f]=∫ab​f(x)dx−k=0∑n​Ak​f(xk​)=...

    1.求积公式余项

    1.1 定义

    R[f]=ab ⁣ ⁣ ⁣f(x)dxk=0nAkf(xk)=Kf(m+1)(η),(1)R[f]=\int_a^b\!\!\!f(x)\mathrm{d}x-\sum_{k=0}^nA_kf(x_k)=Kf^{(m+1)}(\eta ),(1)
    Kf(x),η(a,b).其中K为不依赖f(x)的待定参数,\eta \in (a,b).f(x)m,这个结果表明当f(x)是次数小于等于m的多项式时,
    f(m+1)(x)=0,R[f]=0,ab ⁣f(x)dxk=0nAkf(xk)由于f^{(m+1)}(x)=0,故此时R[f]=0,即求积公式\int_a^b\!f(x)\mathrm{d}x\approx\sum\limits_{k=0}^nA_kf(x_k)
    .f(x)=xm+1,f(m+1)(x)=(m+1)!,(1)Rn(f)0,精确成立.而当f(x)=x^{m+1}时,f^{(m+1)}(x)=(m+1)!,(1)式左端R_n(f)\neq 0,故可求得
    K=1(m+1)![ab ⁣ ⁣ ⁣xm+1dxk=0nAkxm+1]=1(m+1)![1(m+2)(bm+2am+2)k=0nAkxm+1]\begin{aligned} K&=\frac{1}{(m+1)!}\left[\int_a^b\!\!\!x^{m+1}\mathrm{d}x-\sum_{k=0}^nA_kx^{m+1}\right]\\ &=\frac{1}{(m+1)!}\left[\frac{1}{(m+2)}(b^{m+2}-a^{m+2})-\sum_{k=0}^nA_kx^{m+1}\right] \end{aligned}

    1.2 Python实现求积公式余项

    from sympy.abc import x
    from sympy import integrate, exp, diff
    import numpy as np
    from math import factorial
    
    
    def quadrature_reminder(a, b, equal_segment, m, intergrand, a_k: np.ndarray):
        """
        实现[a,b]代数精度为m的求积公式余项表达式
        :param a: 区间左端点即起点a
        :param b: 区间右端点即终点b
        :param equal_segment: 区间等分数n
        :param m: 代数精度m
        :param intergrand: 被积函数f(x)
        :param a_k: 权系数Ak数组
        :return: 求积公式余项表达式
        """
        step_h = (b - a) / equal_segment
        f_x = x ** (m + 1)
        f_x_k_array = np.array([f_x.subs(x, a + k * step_h) for k in range(equal_segment + 1)])
        return 1 / factorial(m + 1) * (1 / (m + 2) * (b ** (m + 2) - a ** (m + 2)) - np.sum(a_k * f_x_k_array)) * diff(
            intergrand, x, m + 1), a, b
    

    2.牛顿-柯特斯公式

    2.1 定义

    [a,b]n,h=ban,xk=a+kh设将积分区间[a,b]分为n等分,步长h=\dfrac{b-a}{n},选取等距节点x_k=a+kh构造出的插值型求积公式
    In=(ba)k=0nCk(n)f(xk),I_n=(b-a)\sum_{k=0}^nC_k^{(n)}f(x_k),
    牛顿-柯特斯(NewtownCotes)公式,Ck(n)柯特斯系数.称为\textbf{牛顿-柯特斯}(\mathrm{Newtown-Cotes})\textbf{公式},式中C_k^{(n)}称为\textbf{柯特斯系数}.
    Ak=ablk(x)dx,k=0,1,,n.x=a+th,根据求积系数A_k=\int_a^bl_k(x)\mathrm{d}x,k=0,1,\ldots,n.引进变换x=a+th,则有
    Ck(n)=hba0nj=0jkntjkjdt=(1)nknk!(nk)!0nj=0jkn(tj)dt.C_k^{(n)}=\frac{h}{b-a}\int_0^n\prod_{j=0\\j\neq k}^n\frac{t-j}{k-j}\mathrm{d}t=\frac{(-1)^{n-k}}{nk!(n-k)!}\int_0^n\prod_{j=0\\j\neq k}^n(t-j)\mathrm{d}t.
    n=2,辛普森(Simpson)公式:特别地,当n=2时,相应的求积公式称为\textbf{辛普森}(\mathrm{Simpson})\textbf{公式}:
    S=ba6[f(a)+4f(a+b2)+f(b)].S=\frac{b-a}{6}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right].

    2.2 Python实现牛顿-柯特斯公式

    def newtown_cotes_integrate(equal_segment, interval_start, interval_end, symbol_t, f_x):
        """
        实现牛顿-柯特斯(Newton-Cotes)插值型积分
        :param symbol_t: 引进变换x=a+t*h后的积分变量
        :param equal_segment: 区间等分数n
        :param interval_start: 区间左端点即起点
        :param interval_end: 区间右端点即终点
        :param f_x: 被积函数 intergrand
        :return:牛顿-柯特斯(Newton-Cotes)插值积分
        """
        step_h = (interval_end - interval_start) / equal_segment
        f_k_array = np.array([f_x.subs(x, interval_start + k * step_h) for k in range(equal_segment + 1)])
        return (interval_end - interval_start) * np.sum(cotes_coefficient(equal_segment, symbol_t) * f_k_array)
    
    
    def cotes_coefficient(equal_segment, symbol_t):
        """
        实现牛顿-柯特斯(Newton-Cotes)插值型求积公式的柯特斯系数
        :param equal_segment: 区间等分数n,其中1<=n<=7,n>7的牛顿-柯特斯公式计算不稳定,不使用.
        :param symbol_t: 引进变换x=a+t*h后的积分变量t
        :return: 柯特斯系数数组
        """
        if equal_segment not in range(1, 8):
            raise ValueError("Cotes coefficient must be an integer between 1 and 7")
        c_k_array = np.zeros(equal_segment + 1, dtype=np.float64)
        for k in range(equal_segment + 1):
            index = list(range(equal_segment + 1))
            index.pop(k)
            intergrand = np.prod(symbol_t - np.array(index))
            integral = integrate(intergrand, (symbol_t, 0, equal_segment))
            c_k_array[k] = (-1) ** (equal_segment - k) / (
                    equal_segment * factorial(k) * factorial(equal_segment - k)) * integral
        return c_k_array
    

    3.复合梯形公式

    3.1 定义

    [a,b]n,xk=a+kh,h=ban,k=0,1,,n,将区间[a,b]划分为n等份,分点x_k=a+kh,h=\dfrac{b-a}{n},k=0,1,\ldots,n,[xk,xk+1](k=0,1,,n1)在每个子区间[x_k,x_{k+1}](k=0,1,\ldots,n-1)上
    ab ⁣f(x)dxba2[f(a)+f(b)],采用梯形公式\int_a^b\!f(x)\mathrm{d}x\approx\dfrac{b-a}{2}[f(a)+f(b)],则得
    Tn=h2k=0n1[f(xk)+f(xk+1)]=h2[f(a)+2k=1n1f(xk)+f(b)],T_n=\frac{h}{2}\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})]=\frac{h}{2}\left[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)\right],
    称为复合梯形公式

    3.2 Python实现复合梯形公式

    def com_trapz(a, b, equal_segment, intergrand):
        """
        实现复合梯形求积公式
        :param a: 区间左端点即起点a
        :param b: 区间右端点即终点b
        :param equal_segment: 区间等分数n
        :param intergrand: 被积函数f(x)
        :return: 复合梯形求积积分
        """
        step_h = (b - a) / equal_segment
        f_x_k_array = np.array([intergrand.subs(x, a + k * step_h) for k in range(1, equal_segment)])
        return step_h / 2 * (intergrand.subs(x, a) + 2 * np.sum(f_x_k_array) + intergrand.subs(x, b))
    

    4.复合辛普森公式

    4.1 定义

    [a,b]n,[xk,xk+1](k=0,1,,n1)将区间[a,b]划分为n等份,在每个子区间[x_k,x_{k+1}](k=0,1,\ldots,n-1)上
    ,xk+1/2=xk+12h,采用辛普森公式,若记x_{k+1/2}=x_k+\frac{1}{2}h,则得
    Sn=h6k=0n1[f(xk)+4f(xk+1/2)+f(xk+1)]=h6[f(a)+4k=0n1f(xk+1/2)+2k=1n1f(xk)+f(b)],\begin{aligned} S_n&=\frac{h}{6}\sum_{k=0}^{n-1}[f(x_k)+4f(x_{k+1/2})+f(x_{k+1})]\\ &=\frac{h}{6}\left[f(a)+4\sum_{k=0}^{n-1}f(x_{k+1/2})+2\sum_{k=1}^{n-1}f(x_k)+f(b)\right], \end{aligned}
    称为复合辛普森公式

    4.2 Python实现复合辛普森公式

    def com_simpson(a, b, equal_segment, intergrand):
        """
        实现复合辛普森求积公式
        :param a: 区间左端点即起点a
        :param b: 区间右端点即终点b
        :param equal_segment: 区间等分数n
        :param intergrand: 被积函数f(x)
        :return: 复合辛普森求积积分
        """
        step_h = (b - a) / equal_segment
        f_x_k1_array = np.array([intergrand.subs(x, a + (k+0.5) * step_h) for k in range(equal_segment)],
                                dtype=np.float64)
        f_x_k2_array = np.array([intergrand.subs(x, a + k * step_h) for k in range(1, equal_segment)],
                                dtype=np.float64)
        return step_h / 6 * (
                intergrand.subs(x, a) + 4 * np.sum(f_x_k1_array) + 2 * np.sum(f_x_k2_array) + intergrand.subs(x, b))
    

    5.测试

    if __name__ == '__main__':
        # 辛普森公式及其余项表达式测试成功,来源详见来源详见李庆扬数值分析第5版P135,e.g.4
        # 0.632333680003663
        inter_grand = exp(-x)
        print("辛普森公式积分为:{}".format(newtown_cotes_integrate(2, 0, 1, x, inter_grand)))
        print("辛普森公式为:{}".format(quadrature_reminder(a=0, b=1, equal_segment=2, m=3, intergrand=exp(-x), a_k=cotes_coefficient(2, x))[0]))
        # 复合梯形公式和复合辛普森公式测试成功,来源详见来源详见李庆扬数值分析第5版P135,e.g.2(1)
        f_x = x / (4 + x ** 2)
        print("复合梯形公式积分:{}".format(com_trapz(a=0, b=1, equal_segment=8, intergrand=f_x)))  # 0.111402354529548
        print("复合辛普森公式积分:{}".format(com_simpson(a=0, b=1, equal_segment=8, intergrand=f_x)))  # 0.111571813252631
    

    6.运行结果

    在这里插入图片描述

    展开全文
  • Python实现复合辛普森求积公式

    千次阅读 2016-01-10 21:44:45
    Python实现复合辛普森求积公式# -*- coding:utf-8 -*- import math def simpr(f, a, b, n): """simpr函数为用复合辛普森公式求积分 f是被积函数 a,b分别为积分的上下限 n是子区间的个数 s是梯形总面积,即所求...
  • 数值积分:梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式1.问题描述微积分方法求积有很大的局限性,当碰到被积函数很复杂时,找不到相应的原函数。积分值在几何上可解释为由 x=a,x=b,y=0和y=...
  • //实验四 复合梯形和辛普森公式 #include<iostream> #include<cmath> using namespace std; class formula{ public: formula(double a,double b,int n); double function(double x); double ...
  • 高斯求积公式的PPT课件。计算方法。 熟悉复合梯形公式、复合抛物线公式及其余项; 熟悉Newton-Cotes公式; 熟悉代数精度法构造求积公式的思想; 熟悉当权为1区间为[-1,1]时的Guass求积公式; 了解变步长梯形公式和...
  • ----------------------个人作业,如果有后辈的作业习题一致,可以参考学习,一起交流,请勿直接copy ··复合抛物线公式: ··龙贝格公式: ...%复合梯形公式 function y=funct...
  • 将积分区间分为若干份, 在每一个“小区间”上用低阶梯形求积公式可得复合梯形公式的收敛阶为2阶。matlab程序function I = ftrapz(fun,a,b,n)%fun,a,b,n分别为被积分函数、积分下限、积分上限、积分区间数目h = (b-a)...
  • 复合梯形公式的提出: 1.首先,什么是梯形公式: 梯形公式表明:f(x)在[a,b]两点之间的积分(面积...但高次的缺陷是当次数大于8次,求积公式就会不稳定。因此,我们用于数值积分的牛顿-科特斯公式通常是一次的梯形...
  • 佛山科学技术学院 实 验 报 告 课程名称 数值分析 实验项目 ... 2学会复合梯形复合Simpson和龙贝格积分公式的编程与应用 3探索二重积分在矩形区域的数值积分方法 二实验要求 按照题目要求完成实验内容 写出相应的M
  • 1、复合梯形公式: #include"stdio.h" #include"math.h" double f1(double x) { return x/(4+x*x); } double f2(double x) { return sqrt(4-(sin(x)*sin(x))); } double Tixing(f,a,b,...
  • 由于其采用的是逐次分半计算,后一次计算是对前一次近似结果的修正,因此相对于辛普森和科特斯求积方法精度更高;并且其前一次分割计算的函数值在分半之后还可以被继续使用,因此提升了计算的效率,是一种精度较高的...
  • 假设被函数为   f x ,积分区间为   , a b ,把区间   , a b 等分成 n 个小区间, 各个区间的长度为 h ,即   / h b a n   ,称之为“步长” 。根据定积分...
  • 在实际应用中,为了既能提高结果的精度,又使算法简便且易在...由此得到的一些具有更大实用价值的数值求积公式统称为复合求积公式。 设f(x)在区间[a,b]上有积分。 则其k步积分递推为:T2k = 1/2T2k-1 + (b-a...
  • 复合梯形积分公式及余项 复合Simpsion积分公式及余项 例题
  • 令Tn为将[a,b]划分n等分的复合梯形求积公式,h =(b-a)/n为小区间的长度。h/2类似于梯形公式中的(b-a)/2注意:这里的k+1是下标通过研究我们发现:T2n与Tn之间存在一些递推关系。注意:这里的k+1/2是下标。并且其中的h...
  • 龙贝格积分

    千次阅读 2013-07-21 17:07:29
    对于不易直接用积分公式计算的原函数,通常用复合梯形求积公式或复合抛物线求积公式等方法,但这些方法精度不高,收敛的速度缓慢。为了提高收敛速度,减少计算量,人们寻求其他方法. Romberg方法也称为逐次分半加速...

空空如也

空空如也

1 2
收藏数 31
精华内容 12
关键字:

复合梯形求积公式