-
梯形求积公式 和 复合梯形求积公式 Matlab 实现
2017-11-04 18:22:52梯形求积公式 和 复合梯形求积公式 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)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)
-
matlab 复合梯形求积公式
2017-03-02 21:41:31matlab 复合梯形求积公式,.M文件,可直接运行出结果。 -
【C语言基础】利用复合梯形求积公式计算定积分
2020-04-03 01:22:38【C语言基础】利用复合梯形求积公式计算定积分 一、复合梯形求积公式 这是数值分析中一种求解定积分的近似方法。适用于被积函数的原函数不能用初等函数表示的情况。 基本思路 将被积函数 f(x)与x轴围成的区域分成n...【C语言基础】利用复合梯形求积公式计算定积分
一、复合梯形求积公式
这是数值分析中一种求解定积分的近似方法。适用于被积函数的原函数不能用初等函数表示的情况。
基本思路
将被积函数 f(x)与x轴围成的区域分成n个梯形,把n个梯形面积求和得到积分的近似值。若精度不满足需要,则可以将每个区间再等分一次,得到2n+1个等分区间,然后再求和,直到精度满足需要。
公式
其中:n 为等分区间数
h=(b-a)/2 (积分步长)
C语言实现算法
以下面积分为例:
首先写出计算被积函数值的函数
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()是用来计算绝对值的函数
-
复合梯形公式和--复合辛普森求积公式----变长梯形公式-------龙贝格方法 ----------matlab
2018-12-06 21:44:56用复合梯形公式和复合辛普森求积公式计算。 复合梯度代码 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function T_n=fht(a,b,n) h=(b-a)/n; for k=0:n x(k+1)=a+k*h; if...- 实验内容
- 计算圆周率π,要求误差小于10-8。
- 用复合梯形公式和复合辛普森求积公式计算。
复合梯度代码 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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
-
数值积分-求积公式余项,牛顿-柯特斯公式,辛普森公式,复合梯形公式,复合辛普森公式
2021-01-06 22:48:581.求积公式余项 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]=∫abf(x)dx−k=0∑nAkf(xk)=...1.求积公式余项
1.1 定义
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 定义
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 定义
称为复合梯形公式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 定义
称为复合辛普森公式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:45Python实现复合辛普森求积公式# -*- coding:utf-8 -*- import math def simpr(f, a, b, n): """simpr函数为用复合辛普森公式求积分 f是被积函数 a,b分别为积分的上下限 n是子区间的个数 s是梯形总面积,即所求... -
数值积分: 梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式
2018-03-09 21:48:49数值积分:梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式1.问题描述微积分方法求积有很大的局限性,当碰到被积函数很复杂时,找不到相应的原函数。积分值在几何上可解释为由 x=a,x=b,y=0和y=... -
c++实现复合求积公式 和辛普森公式
2020-11-29 16:26:37//实验四 复合梯形和辛普森公式 #include<iostream> #include<cmath> using namespace std; class formula{ public: formula(double a,double b,int n); double function(double x); double ... -
高斯求积公式Newton-Cotes公式
2010-01-14 00:14:31高斯求积公式的PPT课件。计算方法。 熟悉复合梯形公式、复合抛物线公式及其余项; 熟悉Newton-Cotes公式; 熟悉代数精度法构造求积公式的思想; 熟悉当权为1区间为[-1,1]时的Guass求积公式; 了解变步长梯形公式和... -
数值计算实验报告---复合求积公式(梯形、抛物线、龙贝格)、导数求值(三点、四点、五点公式)...
2017-04-10 18:44:00----------------------个人作业,如果有后辈的作业习题一致,可以参考学习,一起交流,请勿直接copy ··复合抛物线公式: ··龙贝格公式: ...%复合梯形公式 function y=funct... -
复化梯形公式matlab程序_复合梯形公式(compound trapezoidal formula)
2020-11-24 10:11:41将积分区间分为若干份, 在每一个“小区间”上用低阶梯形求积公式可得复合梯形公式的收敛阶为2阶。matlab程序function I = ftrapz(fun,a,b,n)%fun,a,b,n分别为被积分函数、积分下限、积分上限、积分区间数目h = (b-a)... -
python实现数值分析之龙贝格求积公式
2018-10-24 16:34:23复合梯形公式的提出: 1.首先,什么是梯形公式: 梯形公式表明:f(x)在[a,b]两点之间的积分(面积...但高次的缺陷是当次数大于8次,求积公式就会不稳定。因此,我们用于数值积分的牛顿-科特斯公式通常是一次的梯形... -
MATLAB数值分析实验二(复合梯形、辛普森和龙贝格求积,以及二重积分计算等).doc
2020-12-16 11:18:55佛山科学技术学院 实 验 报 告 课程名称 数值分析 实验项目 ... 2学会复合梯形复合Simpson和龙贝格求积分公式的编程与应用 3探索二重积分在矩形区域的数值积分方法 二实验要求 按照题目要求完成实验内容 写出相应的M -
复合求积法(梯形和Simpson)
2007-03-06 22:04:251、复合梯形公式: #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,... -
复合梯形公式matlab代码_MATLAB龙贝格积分算法
2021-01-27 18:24:15由于其采用的是逐次分半计算,后一次计算是对前一次近似结果的修正,因此相对于辛普森和科特斯求积方法精度更高;并且其前一次分割计算的函数值在分半之后还可以被继续使用,因此提升了计算的效率,是一种精度较高的... -
C语言复合梯形公式实现定积分
2014-03-07 14:51:00假设被积函数为 f x ,积分区间为 , a b ,把区间 , a b 等分成 n 个小区间, 各个区间的长度为 h ,即 / h b a n ,称之为“步长” 。根据定积分... -
复合梯形法的递推算式
2017-11-13 12:20:00在实际应用中,为了既能提高结果的精度,又使算法简便且易在...由此得到的一些具有更大实用价值的数值求积公式统称为复合求积公式。 设f(x)在区间[a,b]上有积分。 则其k步积分递推为:T2k = 1/2T2k-1 + (b-a... -
数值积分22 - 复合求积方法(前面是选有限个点,作n次插值积分。这里是把区间分成n份后,每个区间作n次插值...
2021-01-10 11:54:29复合梯形积分公式及余项 复合Simpsion积分公式及余项 例题 -
python辛普森积分_Python龙贝格法求积分实例
2020-12-02 18:43:28令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方法也称为逐次分半加速...
收藏数
31
精华内容
12
-
ROSE-HA-V8.9+Win2008+SQL2008双机配置详细指南(图文).pdf
-
2021-02-27
-
MMM 集群部署实现 MySQL 高可用和读写分离
-
深究字符编码的奥秘,与乱码说再见
-
libFuzzer视频教程
-
JS//DOM(先占坑)
-
svg_pnf_Selector-源码
-
第二章 分支程序结构设计——作业-答案.html
-
抖音任务点赞平台源码.zip
-
零基础极简以太坊智能合约开发环境搭建并开发部署
-
NFS 实现高可用(DRBD + heartbeat)
-
C++代码规范和Doxygen根据注释自动生成手册
-
【考研打卡】Day027——2020.02.27
-
p != np,这次真的被证明了
-
jupyter的使用.txt
-
Xyplayer X3.9.3正式版.rar
-
JS//异步 单线程(先占坑)
-
qengine:基于查询的处理引擎-源码
-
JS//日期 Math 数组API 对象API(先占坑)
-
RTP协议