第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。
- 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是状态模型文件名称，字符串类型