第3-1章 连续系统数值积分法仿真Matlab编程
数值积分的MATLAB实现
数值积分的MATLAB实现数值积分是通过数值方法计算定积分的近似值。
MATLAB是一种功能强大的数值计算软件，提供了多种函数和工具箱用于数值积分的实现。
在MATLAB中，常用的数值积分方法包括梯形法则、辛普森法则和龙贝格法。
梯形法则是最简单的数值积分方法之一、它的基本思想是将要积分的区间划分成多个小的梯形并计算每个梯形的面积，然后将这些面积相加得到最终的近似积分值。
在MATLAB中，可以使用trapz函数进行梯形法则的计算。
例如，要计算函数sin(x)在区间[0, pi]的积分，可以使用以下代码：```MATLABx = linspace(0, pi, 1000); % 在[0, pi]区间生成1000个等间隔的点y = sin(x); % 计算函数sin(x)在每个点的值integral_value = trapz(x, y) % 使用梯形法则进行数值积分```辛普森法则是一种更精确的数值积分方法，它使用二次多项式来逼近被积函数。
在MATLAB中，可以使用simpson函数进行辛普森法则的计算。
例如，上面例子中的积分可以改用辛普森法则进行计算：```MATLABintegral_value = simpson(x, y) % 使用辛普森法则进行数值积分```龙贝格法是一种高效的自适应数值积分方法，它通过逐步加密网格和逼近函数来提高积分的精度。
在MATLAB中，可以使用quad和quadl函数进行龙贝格法的计算。
例如，计算函数sin(x)在区间[0, pi]的积分：```MATLAB```除了上述方法外，MATLAB还提供了许多其他的数值积分函数和工具箱，用于处理不同类型的积分问题。
例如，int和integral函数可以用于处理多重积分和奇异积分。
Symbolic Math Toolbox中的函数可以用于计算符号积分。
需要注意的是，数值积分是一种近似方法，计算结果的误差与划分区间的精细程度有关。
数值积分的Matlab实现研究
3.指令所采用的算法是将一般的积分区域映射到矩形区域, 然后利用自适应Lobatton算法进行计算.
例题:计算
1
2 1 x 2 2
( x 2 y 2 ) sin(x y 2 )dx
可用如下指令:
quad2d (@(x, y)(x.^2 y.^2).* sin(x y.^2), ...1,2, @( x) x, @( x)1 x.^2)
其中:tic,toc是秒表计时命令,tic表示秒表计时开始,toc表 示秒表计时结束,运行花费时间输出格式为“elapsed_time=”, 单位为秒。 2、被积函数用内嵌函数表示
sym s y 2; y 2 inline('1. /((x 0.3)^2 0.01) 1. /((x 0.9)^2 0.04) 6' );
…*(x.^2+y.^2),x,1+x.^2),x),0,inf)
运行结果:0.6427 例2.计算三重积分 2 2 ( xz y ) ln( 1 x z)dv
其中积分区域由xoy坐标面 与旋转抛物面z=16-x^2-y^2所 围成的立体区域.
4 x 4; 16 x 2 y 16 x 2 ; 0 z 16 x 2 y 2
二、一般数值积分的讨论与研究 由于数值积分非常复杂,如振荡积分,无界区域的积分,无 界函数的积分,高维积分等,目前matlab系统还没有直接提供 一般区域上的三元及三元以上函数积分的指令,这需要人们 从问题出发,利用matlab系统现有的指令,探索求解一般区域 的积分问题. 以具体例题进行讨论: 例1
如 : dz dy ( x 2 y 2 ) ln(x y z )dx
第3章 连续系统数字仿真的基本算法
(2.18)
h y k 1 y k (k1 3k3 ) 4 k1 f (t k , y k ) h h k 2 f (t k , y k k1 ) 3 3 k3 f (t k 2 h, y k 2 hk2 ) 3 3
上一页
下一页
返回
表2.1 f 的计算次数与算法精度阶数的关系
每步计算f 的次 数 算法精度阶数 2 3 4 4 5 6 r-2
2
3
4
5
6
7
r≥8
由此可见，RK4法有其优越性。
上一页
下一页
返回
3.1.4 微分方程数值积分的矩阵分析
对于一阶向量微分方程及初值问题
y f (t , y) y( t0 ) y0
(2.16)
W1 W2 1 1 W2 c2 2 1 W2 a21 2
取c2=1，有
(2.17)
W1 W2
上一页
1 2
a21 1
下一页 返回
从而，有RK2法
h y k 1 y k 2 (k1 k 2 ) k1 f (t k , y k ) k 2 f (t k h, y k hk1 )
第3章 连续系统数字仿真的基 本算法
2.1 数值积分算法 2.2 数值积分算法的基本分析 2.3 连续系统仿真的离散相似算法 2.4 常用快速数字仿真算法 2.5 实时数字仿真算法 小结
3.1 数值积分算法
3.1.1 3.1.2 3.1.3 3.1.4 数值积分算法的基本原理 欧拉法 龙格-库塔法 微分方程数值积分的矩阵分析
上一页 下一页 返回
数值积分的matlab实现
实验10 数值积分实验目的：1．了解数值积分的基本原理； 2．熟练掌握数值积分的MATLAB 实现； 3．会用数值积分方法解决一些实际问题。
实验内容：积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。
同微分一样，在《微积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。
此外有些函数虽然有解析式，但其原函数不是初等函数，所以仍然得不到积分的精确值，如不定积分⎰10 d sin x x x。
这时我们一般考虑用数值方法计算其近似值，称为数值积分。
10.1 数值微分简介设函数()y f x =在*x 可导，则其导数为hx f h x f x f h )()(lim )(**0*-+='→ （10.1）如果函数()y f x =以列表形式给出（见表10-1），则其精确值无法求得，但可由下式求得其近似值hx f h x f x f )()()(***-+≈' （10.2）表 10-1一般的，步长h 越小，所得结果越精确。
（10.2）式右端项的分子称为函数()y f x =在*x 的差分，分母称为自变量在*x 的差分，所以右端项又称为差商。
数值微分即用差商近似代替微商。
常用的差商公式为：000()()()2f x h f x h f x h +--'≈（10.3）hy y y x f 243)(2100-+-≈' （10.4）hy y y x f nn n n 234)(12+-≈'-- （10.5）其误差均为2()O h ，称为统称三点公式。
10.2 数值微分的MATLAB 实现MATLAB 提供了一个指令求解一阶向前差分，其使用格式为： dx=diff(x) 其中x 是n 维数组，dx 为1n -维数组[]21321,,,n x x x x x x ---,这样基于两点的数值导数可通过指令diff(x)/h 实现。
matlab 连续积分和非连续微分
Matlab是一种功能强大的数学软件，广泛应用于工程、科学和数学领域。
其中，连续积分和非连续微分是其重要功能之一。
本文将对Matlab中的连续积分和非连续微分进行详细介绍。
1. 连续积分在Matlab中，连续积分是指对一个连续函数进行积分运算。
积分是微积分的一个重要概念，用于求函数曲线下的面积以及确定函数的反导数。
在Matlab中，可以使用int函数进行连续积分运算。
对于函数f(x) = x^2，我们可以使用Matlab进行连续积分运算，其代码如下：```matlabsyms xF = int(x^2, x)```上述代码中，我们首先使用syms声明x为符号变量，然后使用int函数对x^2进行积分运算，得到的结果为F = 1/3*x^3。
Matlab还提供了更加丰富的连续积分功能，包括定积分、不定积分、多重积分等。
通过这些功能，用户可以方便地进行各种积分运算，满足不同的数学求解需求。
2. 非连续微分在实际问题中，有些函数是不连续的，这时需要对这些不连续函数进行微分运算。
在Matlab中，可以使用diff函数对不连续函数进行微分运算。
对于函数f(x) = |x|，我们可以使用Matlab进行微分运算，其代码如下：```matlabsyms xf = abs(x)df = diff(f, x)```上述代码中，我们首先使用syms声明x为符号变量，然后使用diff 函数对| x |进行微分运算，得到的结果为df = sign(x)。
除了对不连续函数的微分运算外，Matlab还提供了更多的微分功能，包括高阶微分、偏微分、方程微分等。
这些功能可以帮助用户轻松解决各种微分运算问题。
总结：Matlab作为一种强大的数学软件，提供了丰富的积分和微分功能，包括连续积分、非连续微分等。
利用这些功能，用户可以高效地进行各种数学求解和模拟，满足实际问题的需求。
希望本文对Matlab的积分和微分功能有所帮助，同时也希望读者能够更深入地了解和应用这些功能，进一步提升数学建模和求解的能力。
数值分析MATLAB编程——数值积分法
数值分析MATLAB编程——数值积分法1、调用函数--f.Mfunction y=f(x)%------------------------------------------------------------函数1 y=sqrt(4-sin(x)*sin(x));%------------------------------------------------------------函数2 %y=sin(x)/x;%if x==0% y=0;%end%------------------------------------------------------------函数3 %y=exp(x)/(4+x*x);%------------------------------------------------------------函数4 %y=(log(1+x))/(1+x*x);2、复合梯形公式--tixing.M%复合梯形公式clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');h=(b-a)/n;x=a:h:b;T=0;for k=1:n;T=0.5*h*(f(x(k))+f(x(k+1)))+T;endT=vpa(T,8)3、复合Simpson公式--simpson.M%复合Simpson公式clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');h=(b-a)/n;x=a:h:b;S=0;for k=1:n;xx=(x(k)+x(k+1))/2;S=(1/6)*h*(f(x(k))+4*f(xx)+f(x(k+1)))+S;endS=vpa(S,8)4、Romberg算法--romberg.M%Romberg算法clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');num=0:n;R=[num'];h=b-a;T=h*(f(a)+f(b))/2;t(1)=T;for i=2:n+1;u=h/2;H=0;x=a+u;while x<b;H=H+f(x);x=x+h;endt(i)=(T+h*H)/2;T=t(i);h=u;endR=[R,t'];for i=2:n+1for j=n+1:-1:1if j>=it(j)=(4^(i-1)*t(j)-t(j-1))/(4^(i-1)-1);elset(j)=0;endendR=[R,t'];endR=vpa(R,8)R(n,n)5、变步长算法（以复化梯形公式为例）--tixing2.M%复合梯形公式，确定最佳步长format longclear alla=input('请输入积分下限:');b=input('请输入积分上限:');eps=input('请输入误差：');k=1;T1=(b-a)*(f(a)+f(b))/2;T2=(T1+(b-a)*(f((a+b)/2)))/2; while abs((T1-T2)/3)>=epsM=0;n=2^k;h=(b-a)/n;T1=T2;x=a:h:b;for i=1:n;xx=(x(i)+x(i+1))/2;M=M+f(xx);endT2=(T1+h*M)/2;k=k+1;endT=vpa(T2,8)n=2^k。
3-1连续系统数值积分仿真方法学
21
龙格-库塔法的特点
A 单步法
可自启动，即知道初值后，直接从微分方程
初值开始计算下去。 B 不同步骤的步长可以不一样，但每一步需 同一步长来计算系数。 C 计算YK+1分两部分：求YK；步长乘各系数 的加权平均值
2019/2/6 22
D 精度取决于h的大小及解决方法 当精度相同时，四阶的计算量是二阶的两倍，
2019/2/6
k3 f yk k 2 h / 2, t k h / 2
20
k1 f yk , t k
yk 1 yk hk1 3k 2 3k3 k 4 / 8
k 2 f yk k1h / 3, t k h / 3
k3 f yk k 2 h k1h / 3, t k 2h / 3 k 4 f yk k1h k 2 h k3 h, t k h
B—n*1控制矩阵;C—1*n输出矩阵
U(t)—输入变量;y(t)—输出变量
采用四阶龙格库塔法进行仿真计算
2019/2/6
33
龙格库塔公式
xk 1
k1 f yk , t k
xk hk1 2k 2 2k3 k 4 / 6
k 2 f yk k1h / 2, t k h / 2 k 4 f yk k3 h, t k h
2019/2/6 30
二、面向方程的系统仿真
1 状态方程的仿真程序
2 传递函数
3 微分方程
2019/2/6
31
1 状态方程的仿真程序
设系统的状态方程和输出方程：
t Axt But x yt Cxt
数值积分算法与MATLAB实现陈悦5133201讲解
东北大学秦皇岛分校数值计算课程设计报告数值积分算法及MATLA实现学院数学与统计学院专业信息与计算科学学号5133201姓名_____________ 陈世_____________指导教师姜玉山张建波成绩教师评语:指导教师签字:2015年07月14日1 绪论数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现•而数值分析主要研究数值计算•现科学技术的发展与进步提出了越来越多的复杂的数值计算问题，这些问题的圆满解决已远人工手算所能胜任，必须依靠电子计算机快速准确的数据处理能力•这种用计算机处理数值问题的方法，成为科学计算•今天，科学计算的应用范围非常广泛，天气预报、工程设计、流体计算、经济规划和预测以及国防尖端的一些科研项目，如核武器的研制、导弹和火箭的发射等，始终是科学计算最为活跃的领域•1.1数值积分介绍数值积分是数值分析的重要环节，实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的•另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解•由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题•对微积分学做出杰出贡献的数学大师，如I.牛顿、L.欧拉、C.F.高斯、拉格朗日等人都在数值积分这个领域作出了各自的贡献，并奠定了这个分支的理论基础•构造数值积分公式最通常的方法是用积分区间上的n次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式•特别在节点分布等距的情形称为牛顿-科特斯公式，例如梯形公式(Trapezoidal Approximations)与抛物线公式(Approximatens Using Parabolas)就是最基本的近似公式•但它们的精度较差•龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式(Rhomberg Integration).当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分•数值积分还是微分方程数值解法的重要依据•许多重要公式都可以用数值积分方程导出•现探讨数值积分算法以及运用MATLAB软件的具体实现1.2 MATLAB 软件MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分.MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）.是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境.它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平•MATLAB和Mathematics Maple并称为三大数学软件.它在数学类科技应用软件中在数值计算方面首屈一指.MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域.MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件.在新的版本中也加入了对C，FORTRAN, C++，JAVA的支持.2数值积分的基本概念一般的，我们可以在区间la,b 1上适当选取某些节点兀，然后用f X k的加权平均得到平均高度f 的近似值，这样构造出的求积公式具有下列形式：a n八A k f X k ，式中X k称为求积节点；A k称为求积系数，亦称伴随节点X k的k =0b f x dx ：权.权A k仅仅与节点X k的选取有关，而不依赖于被积函数 f x的具体形式.2.1代数精度的概念如果某个求积公式对于次数不超过m的多项式均能准确的成立，但对于m T次多项式就不准确成立，则称该求积公式具有m次代数精度（或代数精确值）一般地，欲使求积公式具有m次代数精度，只要令它对于 f x =1,x,|||,x m都能准确成立，这就要求：乞A k =b-a,1 2 2为A k x k=2(b -a)‘++-- - m 1 t. m+1 m+1 ,送A k x k =一 b -a .、一m+12.2求积公式的余项b n令求积公式的余项为R【f】，其中R i f ]二J a f (x)dx-送A k f(x J;区间【a,b】可以是k=0有限的或无限的•构造求积公式的问题就是确定X j和A j使得R〔f 1在某种意义下尽可能地小.3数值积分方法及MATLAB实现3.1复合辛普森公式3.1.1插值型求积公式用插值多项式L n(X)替换积分广=[f(XpX中的被积函数f(X)，然后计算bI 「L n XdX作为积分的近似值，这样建立的求积公式称为插值型求积公式•用插值多an n项式L n(x)的表达式L n(X)八I k X f X k，代入得I二^ A k f X k，其中：k=0 k=0A = f (X—X。
计算机仿真技术课件3数值积分法在系统仿真中的应用
第三章 数值积分法在系统仿真中的应用
龙格-库塔法
将式(3.3)在tn点展开 泰勒级数
h2 (t n ) (t n ) o(h 3 ) y (t n h) y (t n ) hy y 2
龙格-库塔（RK）法的一般形式为
r
（3-9）
yn1 yn h wi ki
第三章 数值积分法在系统仿真中的应用
为了适用于实时仿真计算，一般经常采用以下方法：
（1）选择Adams多步法。 （2）合理地选择龙格-库塔法计算公式中的系数，使之适用于 实时仿真。
在方程(3 18)式中，令W1 0可得W2 1, c2 1 , 此时， 2 (3 18)式化为
yn 1 yn hk2 k1 f (t n , yn ) k 2 f (t n h , yn h k1) 2 2
（3-38）
F为函数，外部输入为u(t) 。
第三章 数值积分法在系统仿真中的应用
采入un1
采入u (tn h ) 2
计算yn1 , 并输入
计算k 2
计算下 一个k1
tn
tn h 2
t n 1
图3.6
tn1 h 2
tn1 h
RK-2 的计算流程
由于当tn h tn 1时极端即才能真正采入 外部输入un1 , 才具备 计算k2的条件。所以 yn1要到tn1 h 时才能计算出来，并输 入到外 2 部设备，也就是说，计 算机输出要迟后半个计 算步距。
p yn 1 yn h f ( yn , t n ) （3-8） p p h y y [ f ( y , t ) f ( y n n n 1 , t n 1 )] n 1 n 2
连续系统仿真的方法.
第3章 连续系统仿真的方法3.1 数值积分法连续系统数值积分法，就是利用数值积分方法对广微分方程建立离散化形式的数学模型——差分方程，并求其数值解。
可以想象在数学计算机上构造若干个数字积分器，利用这些数字积分器进行积分运算。
在数字计算机上构造数字积分器的方法就是数值积分法，因而数字机的硬件特点决定了这种积分运算必须是离散和串行的。
把被仿真系统表示成一阶微分方程组或状态方程的形式。
一阶向量微分方程及初值为()(),00t Y Y t Y ⎧⎫⎪⎪⎨⎬⎪⎪⎩⎭Y =F =（3-1）其中，Y 为n 维状态向量，F （t ，Y ）为n 维向量函数。
设方程（3-1）在011,,,,n n t t t t t +=…处的形式上的连续解为()()()()n+1n+1t t n+10t t t =Y t +,(),n Y F t Y dt Y t F t Y dt=+⎰⎰（3-2）设 n =()n Y Y t ，令1n n n Y Y Q +=+（3-3）则有：()1n+1t n Y Y +=也就是说，1(,)n nt n t Q F t Y dt +≈⎰（3-4）如果n Y 准确解()n Y t 为近似值，n Q 是准确积分值的近似值，则式（3-4）就是式（3-2）的近似公式。
换句话说，连续系统的数值解就转化为相邻两个时间点上的数值积分问题。
因此，所谓数值解法，就是寻求初值问题（3-1）的真解在一系列离散点12n t t t <…<…上的近似解12,,,n Y Y Y ……，相邻两个时间离散点的间隔1n n n t t +=-h ，称为计算步距或步长，通常取n =h h 为定值。
可见，数值积分法的主要问题归结为对函数(,)F t y 的数值积分问题，即如何求出该函数定积分的近似解。
为此，首先要把连续变量问题用数值积分方法转化成离散的差分方程的初值问题，然后根据已知的初值条件0y ，逐步地递推计算后续时刻的数值解(1,2,)i y i =…。
连续系统数字仿真数值积分法
main() { float x11,x21,x31,x10,x20,x30,R; int i; x10=0;x20=0;x30=0;outputy[0][0]=0; ST=10;DT=0.005; LP=ST/DT;VN=1;R=20; for (i=1;i<=LP;i++) { x11=x10-DT*x30+DT*R; x21=DT/0.5*x10+(1-DT/0.5)*x20; x31=DT/0.1*x20+(1-DT/0.1)*x30; outputy[0][i]=x31; x10=x11;x20=x21;x30=x31; } dispcurve(); } 程序为llx56.c
4．2 梯形法 为了提高仿真精度，离散-再现环节采取图4．3 的形式。 (4-11)
X (k 1) X (k )
T [e(k ) e(k 1)] 2
(4-12)
式(4-12)称为梯形公式，其几何解释如图4．4所示。
[例4．1] 已知一多变量系统的结构框图如图4．5所示，请 用梯形公式对此系统进行仿真，并输出y1、y2的仿真结果。
假设线性系统的状态空间描述为
AX BU X Y CX DU
式中：X为n×1维状态向量；U为r×l维输入 向量；A为n×n维系统矩阵；B为n×r维输入 矩阵；Y为m×1维输出向量；C为m×n维输出 矩阵；D为m×r维传递矩阵。
(4-3) (4-4) (4-5) (4-6)
4.3 龙格-库塔(Runge-Kutta)法
在非实时仿真中，有时需要更高的精度。 4．3节和4．4节中将再介绍两种更精确的方 法及其离散-再现环节。
MATLAB及控制系统仿真数值积分
x(t0 ) f (t , x, u (t ))dt f (t , x, u (t )) dt
x(tn1 ) x(tn ) f (t , x, u (t ))dt
tn t n1
t x n1 xn Qn 希望找到一个近似公式 x( t n1 ) x( t n ) f ( t , x , u( t ))dt t 来表示方程的近似解： t
t n 1
t n 1
tn
f (t , x ,u(t ))dt
tn
1 f (t , x,u(t ))dt h f t n , x n,un f t n 1, x n 1,un 1 2
h 写成差分方程为 x n 1 x n fn fn 1 2
梯形近似及其误差
前向欧拉yn yn1 hf (tn1 , yn1 )
后向欧拉yn yn1 hf (tn , yn )
前向欧拉yn yn1 hf (tn1 , yn1 )
后向欧拉yn yn1 hf (tn , yn )
f(t) 误差 fn fn+1
f(t) fn
误差
fn+1
f(t) 误差 fn fn+1
梯形近似
tn
tn+1
t
f(t) 误差 fn
fn+1
f(t) fn
误差
fn+1
左矩形近似
右矩形近似
tn
tn+1
t
tn
tn+1
t
f(t) 误差 fn
fn+1
梯形近似
tn
tn+1
第三章 matlab的simulink建模与仿真
nonlinear control
4、提供仿真库的扩充和定制功能
5、应用领域
通信与卫星系统 航空航天
生物系统
汽车系统
船舶系统
金融系统
3、simulink在matlab家族中的位置
Stateflow Blockset Toolboxes coder RTW compiler
simulink MATLAB
第三章 matlab的simulink建模 与仿真
3.1 绪论
一、系统与模型
1、系统
系统是指具有某些特定功能，相互联系、相互作 用的元素集合。 系统的两个基本特征：整体性、相关性
对系统的研究从以下三个方面入手：
1）实体：组成系统的元素，对象
2）属性：实体的特征
3）活动：系统状态变化的过程
系统仿真是研究系统的一种重要手段，而系统模 型是仿真所研究的直接对象。 2、系统模型 实体模型：根据相似性建立 模型 数学模型:原始系统数学模型；仿真系统数学模型
连续系统的输入输出方程为： y(t ) u(t ) sin u(t ) u(t)与y(n)的数学关系为： u(t ) y(n), nTs t (n 1)Ts 整个系统的方程描述：
y (t ) u (n) n / 2, n 1,2,3... y (n) u (n) 1, y (n) sin( y (n)),n t n 1
Function&Tables（函数与表库）
表数据选择器（从表中选择数据） 求取输入信号的数学函数值 对输入信号进行内插运算 输入信号的一维线性内插 输入信号的二维线性内插
输入信号的n维线性内插
M函数（对输入进行运算输出结果） 多项式求值 查找输入信号所在范围 S-函数模块 S-函数生成器
计 算 机 仿 真 技 术第三章 数值积分法在系统仿真中的应用
取 c2 1，得
W1=
W2=
1 2
a21 1
故得二阶龙格—库塔法计算公式
yn1
yn
1 2
K1
K2
K1 hf tn , yn K2 hf tn h, yn K1
(3.1.19)
由于(3.1.13)式中只取了 h ，h 2两项，而将h 2 以上的高阶项 忽略了，所以这种计算方法的截断误差正比于h 3。
ytn
h
ytn
hf
tn ,
yn
h h2 2!
f
tn ,
yn
将(3.1.4)式在
Rn
1 2!
h2
f
tn
,
yn
以后截断，即
得(3.1.3)式的欧拉公式，Rn称为局部截断误差，
它与 h2成正比，即 Rn O h2
(3.1.5)
另外，解以t 0开始继续到t tn ，所积Βιβλιοθήκη 的误3.1.2龙格—库塔法
欧拉法是将 y f t, y, yt1 y0在 tn点附近 的 ytn h经台劳级数展开并截去 h2以后各项得到
的一阶一步法，所以精度较低。如果将展开式 (3.1.4)式多取几项以后截断，就得到精度较高的 高阶数值解，但直接使用台劳展开式要计算函数 的高阶导数。龙格—库塔法是采用间接利用台劳 展开式的思路，即用在n个点上的函数值f的线性 组合来代替f的导数，然后按台劳展开式确定其中 的系数，以提高算法的阶数。这样既能避免计算 函数的导数，同时又保证了计算精度。由于龙 格—库塔法具有许多优点，故在许多仿真程序包 中，它是一个最基本的算法之一。
第3-1章 连续系统数值积分法仿真Matlab编程
6
3.2 算法稳定性分析仿真编程
针对下面的系统，求用解析法、欧拉法和RK4分别求解，计算欧 拉法最大允许步长，将步长从0.1逐渐增大，比较三种解的效果。
4 x x x(0) 1
解：1）该系统是稳定的，解析解为
(3.2-1)
x(t ) e4t
2）用欧拉法计算本例时，其步长应该满足
% 输出参数：NewX是这一步计算的新的状态向量
单步数值积分函数只是对微分方程组StateModel进行一步的计算， 计算法由InteMethod参数指定，可以上欧拉法，RK2或RK4。
4
function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel); if InteMethod == 'RK4' k1=eval([StateModel,'(curt,curx,curu)']); k2=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k1,curu)']); k3=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k2,curu)']); k4=eval([StateModel,'(curt+h,curx+h*k3,curu)']); NewX=curx+h*(k1+2*k2+2*k3+k4)/6; elseif InteMethod == 'RK2' k1=eval([StateModel,'(curt,curx,curu)']); k2=eval([StateModel,'(curt+h,curx+h*k1,curu)']); NewX=curx+0.5*h*(k1+k2); else %欧拉法EUL Qk=eval([StateModel,'(curt,curx,curu)']);%eval用来执行一个函数，传递函数 名和函数的输入参数 NewX=curx+h*Qk; end
数值积分的matlab实现
然后编写如下程序,并保存为t2.m
%t2.m
clear;
a=0;b=2*pi;
T=chen_trap(@f2,a,b);
S=chen_simpson(@f2,a,b);
FT8=comtrap(@f2,a,b,8);
FS4=comsimpson(@f2,a,b,8);
quad=Romberg(@f2,a,b,10e-6);
(4)
(5)
1、依照实验原理编写各种公式的的程序.
2、首先在电脑上安装matlab，然后启动matlab，分别建立不同公式的M文件，
实验程序如下；
梯形公式程序代码
程序代码说明
function T=chen_trap
(fun,a,b)；
T=(b-a)*(feval(fun,a)+feval(fun,b))/2；
disp('与精确值比较的误差为:'),
wuchi=[double(abs(Sj-T)),double(abs(Sj-S)),double(abs(Sj-FT8)),double(abs(Sj-FS4)),double(abs(Sj-quad))]
五、实验数据及结果：
运行t1得
t1
积分准确值为:
Sj = 0.492712
disp('与精确值比较的误差为:'),
wuchi=[double(abs(Sj-T)),double(abs(Sj-S)),double(abs(Sj-FT8)),double(abs(Sj-FS4)),double(abs(Sj-quad))]
;
首先编写被积函数的M文件f2.m
function y=f2(x)
(2021年整理)数值积分用matlab实现
数值积分用matlab实现编辑整理：尊敬的读者朋友们：这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（数值积分用matlab实现）的内容能够给您的工作和学习带来便利。
同时也真诚的希望收到您的建议和反馈，这将是我们进步的源泉，前进的动力。
本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快业绩进步，以下为数值积分用matlab实现的全部内容。
东北大学秦皇岛分校数值计算课程设计报告数值积分及Matlab实现学院数学与统计学院专业信息与计算科学学号5133117姓名楚文玉指导教师张建波姜玉山成绩教师评语：指导教师签字：2015年07月14日1 绪论在科研计算中，经常会碰到一些很难用公式定理直接求出精确解的积分问题，对于这类问题，我们一般转化为数值积分问题，用计算机来实现求解问题． 1.1 课题的背景对于定积分()ba f x dx ⎰在求某函数的定积分时,在一定条件下，虽然有牛顿—莱布里茨公式()()()ba I f x dx Fb F a ==-⎰可以计算定积分的值，但在很多情况下的原函数()f x 不易求出或非常复杂．被积函数的原函数很难用初等函数表达出来,例如2sin (),x x f x e x-=等;有的函数()f x 的原函数()F x 存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式．因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的．另外,许多实际问题中的被积函数()f x 往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值．因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算．而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值．微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节数值积分是数学上重要的课题之一，是数值分析中重要的内容之一．随着计算机的出现，近几十年来,对于数值积分问题的研究已经成为一个很活跃的研究领域．现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义．国内外众多学者在数值积分应用领域也提出了许多新方法．在很多实际应用中,只能知道积分函数在某些特定点的取值,比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等．通过这个课题的研究，我们将会更好地掌握运用数值积分算法求出特殊积分函数的定积分的一些基本方法、理论基础；并且通过Matlab 软件编程的实现，应用于实际生活中． 1.2 课题的主要内容框架1。
数值积分算法与MATLAB实现 论
编号：审定成绩：毕业设计（论文）设计（论文）题目：数值积分算法与MATLAB实现学院名称：数理学院学生姓名：专业：数学与应用数学班级：学号：指导教师：答辩组负责人：填表时间：年月摘要在求一些函数的定积分时，由于原函数十分复杂难以求出或用初等函数表达，导致积分很难精确求出，只能设法求其近似值，因此能够直接借助牛顿-莱布尼兹公式计算定积分的情形是不多的。
数值积分就是解决此类问题的一种行之有效的方法。
积分的数值计算是数值分析的一个重要分支；因此，探讨近似计算的数值积分方法是有着明显的实际意义的。
本文从数值积分问题的产生出发，详细介绍了一些数值积分的重要方法。
本文较详细地介绍了牛顿-科特斯求积公式，以及为了提高积分计算精度的高精度数值积分公式，即龙贝格求积公式和高斯-勒让德求积公式。
除了研究这些数值积分算法的理论外，本文还将这些数值积分算法在计算机上通过MATLAB软件编程实现，并通过实例用各种求积公式进行运算，分析比较了各种求积公式的计算误差。
【关键词】数值积分牛顿-科特斯求积公式高精度求积公式 MATLAB软件ABSTRACTWhen the solution of the definite integral of some function values，because the original function is very complex and difficult to find the elementary function expression, the integral is difficult to accurately calculate, only managed to find the approximate value, and the case is small that allows to direct interface with the Newton - Leibniz formula to calculate the definite integral. Numerical integration is an effective method to solve such problems. The numerical integration is an important branch of numerical analysis; therefore, exploring the approximate calculation of the numerical integration method has obvious practical significance. This article departure from the numerical integration problem, described in detail some important numerical integration methods.This paper has introduced detail the Newton - Coates quadrature formula, and in order to improve the calculation accuracy of numerical integration formulas, More precise formulas have Romberg quadrature formulas and the Gauss - Legendre quadrature formula. In addition to the study of these numerical integration algorithm theory, the article also involve what these numerical integration algorithm be programmed by matlab software on the computer, and an example is calculated with a variety of quadrature formulas, finally analysis and comparison to various quadrature formulas calculation error.【Key words】Numerical integration Newton-Cotes quadrature formulaHigh-precisionquadrature formula Matlab software目录前言 (1)第一章牛顿-科特斯求积公式 (2)第一节数值求积公式的构造 (2)第二节复化求积公式 (9)第三节本章小结 (12)第二章高精度数值积分算法 (13)第一节梯形法的递推 (13)第二节龙贝格求积公式 (14)第三节高斯求积公式 (17)第四节高斯-勒让德求积公式 (19)第五节复化两点高斯-勒让德求积公式 (22)第六节本章小结 (23)第三章各种求积公式的MATLAB编程实现与应用 (24)第一节几个低次牛顿-科特斯求积公式的MATLAB实现 (24)第二节复化求积公式的MATLAB实现 (28)第三节龙贝格求积公式的MATLAB实现 (33)第三节高斯-勒让德求积公式的MATLAB实现 (34)第五节各种求积算法的分析比较 (36)第六节本章小结 (38)结论 (39)致谢 (40)参考文献 (41)附录 (43)一、英文原文 (43)二、英文翻译 (53)前 言对于定积分()ba f x dx ⎰，在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式()()()ba I f x dx Fb F a ==-⎰可以计算定积分的值，但在很多情况下()f x 的原函数不易求出或非常复杂。
数值积分法的matlab程序
4、实验步骤：
（1）编写 matlab 程序，每种方法选取 0.01、0.1、0.5*h 临界、h 临界四种步长求数值 解，分别与解析解比较。画出 matlab 程序的流程图。
（2）结果输出：生成数据文件，将四种数值方法的结果分别与解析法的计算结果存在 下 面的 ASCII 码 数据 文件 中，分别 命名为： eulerdata.dat ， eulerpdata.dat ， rkdata.dat ， adamsdata.dat。 已生成，例如下图为龙哥库塔法、h=0.01 时的部分数据文件，从数据可知，该方法的 仿真结果较好，与数值解的差别不大。
k_1=a*y1(:,i)+b*u; k_2=a*y1(:,i-1)+b*u; k_3=a*y1(:,i-2)+b*u; k_4=a*y1(:,i-3)+b*u; x1=x+h/24*(55*k_1-59*k_2+37*k_3-9*k_4); k_5=a*x1+b*u; x=x+h/24*(9*k_5+19*k_1-5*k_2+k_3); y=[y;x']; y1=y'; end t1=0:0.01:10; y2=1.84-4.95*t1.*exp(-1.88*t1)-1.5*exp(-1.88*t1)-0.34*exp(-6.24*t 1); subplot(2,2,k) plot(t,y) hold on plot(t1,y2,'k-.','linewidth',2) xlabel('time(sec)') ylabel('·ùÖµy') legend('y','doty','ddoty','y-shuzhi') if (k==1) title('admas h=0.01') axis([0,10,-1,3]) elseif (k==2) title('admas h=0.1') axis([0,10,-1,3]) elseif (k==3) title('admas h=0.105 axis([0,10,-1,3]) else (k==4) title('admas h=0.21 临界') axis([0,20,-3,3]) end end 0.5*临界 0 0 0 1 ，b 设： a 0 0 。 22.06 27 10 40.6
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
函数w_DigiInteSimu和w_StepIntegral构造了一个数值积分法 仿真的框架，并不涉及具体的系统。
具体的系统由StateModel,ControlFile,OutputFile参个参数决定， 实际上就是三个函数文件名，这三个函数输入输出参数必须遵循特定 的格式，在准备好由这3个函数描述的系统后，调用w_DigiInteSimu 即可进行仿真。 还需要一个调用w_DigiInteSimu进行仿真的程序，它指定模型文 件，指定初始参数，并且对仿真结果绘图。
2
function [t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile, ControlFile) t=[tstart:h:tstop];%t数一个行序列 cntt=size(t,2);%返回列数 y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果 y0=eval([OutputFile,'(tstart,x0,u0)']);%计算初始输出 y(:,1)=y0’;%将cury作为输出的第1列 curx=x0; %当前一步的x curu=u0; %当前一步的u cury=y0; %当前一步的y for i=1:1:cntt-1 curu=eval([ControlFile,'(t(i),h,curx,curu,cury)']);%计算控制时传递的参数：当前 时间，步长，当前状态和输出 curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%单步积分运算 cury=eval([OutputFile,'(t(i),curx,curu)']);%计算输出 y(:,i+1)=cury‘;%将输出加入到输出序列里 end
(3.1-1)
对这样的系统进行仿真，实际上涉及到控制的计算、状态的数值 积分计算和输出的计算。3个函数g,f,h确定后，就可以完整地描 述一个系统。
1
(1)数值积分法仿真框架函数
function [t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile, ControlFile) % 函数功能：用数值积分法仿真 % 输入参数：tstart, tstop,h 分别是起始时间、结束时间和仿真步长，是标量 % x0,u0是状态、输入的初值，都是列向量 % cnty是输出变量的个数 % InteMethod时数值积分方法，可以是'EUL' ,'RK2','RK4' % StateModel是状态模型的文件名 % ControlFile是控制作用的文件名 % OutputFile是系统输出的文件名 % 输出参数：t是仿真结果的时间序列 % y是仿真结果系统的输出序列
1 4 h 1 h 0 . 5
3）RK4步长条件式是一个高阶不等式，无法直接求解，只能用 试探法确定RK4的步长上限。
% 输出参数：NewX是这一步计算的新的状态向量
单步数值积分函数只是对微分方程组StateModel进行一步的计算， 计算法由InteMethod参数指定，可以上欧拉法，RK2或RK4。
4
function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel); if InteMethod == 'RK4' k1=eval([StateModel,'(curt,curx,curu)']); k2=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k1,curu)']); k3=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k2,curu)']); k4=eval([StateModel,'(curt+h,curx+h*k3,curu)']); NewX=curx+h*(k1+2*k2+2*k3+k4)/6; elseif InteMethod == 'RK2' k1=eval([StateModel,'(curt,curx,curu)']); k2=eval([StateModel,'(curt+h,curx+h*k1,curu)']); NewX=curx+0.5*h*(k1+k2); else %欧拉法EUL Qk=eval([StateModel,'(curt,curx,curu)']);%eval用来执行一个函数，传递函数 名和函数的输入参数 NewX=curx+h*Qk; end
6
3.2 算法稳定性分析仿真编程
针对下面的系统，求用解析法、欧拉法和RK4分别求解，计算欧 拉法最大允许步长，将步长从0 x x x (0) 1
解：1）该系统是稳定的，解析解为
(3.2-1)
x(t) e4t
2）用欧拉法计算本例时，其步长应该满足
第三章 Matlab编程(数值积分法仿真)
3.1 连续系统数值积分法仿真编程思路
目的：针对下面的系统编写通用的数值积分法仿真程序
u g (t , x , u , y ) f (t , x , u ) x y h (t , x , u )
给定初值：u0，y0,系统中的向量都采用列向量表示
3
(2)单步数值积分法函数
function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel) %函数功能：单步积分运算,与模型方程有关 % 输入参数：curt是当前时间，h是数值积分步长 % % % curx,curu分别是当前的状态和控制向量 InteMethod是积分算法，字符串类型，可以取'EUL','RK2','RK4' StateModel是状态模型文件名称，字符串类型