数值积分法的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实现romberg积分法
一、概述在数值分析中，求解定积分是一项重要的任务。
传统的数值积分方法包括梯形法则、辛普森法则等。
而Romberg积分法，是一种更加精确的数值积分方法，它通过不断增加区间的细分，逐步提高数值积分的精度。
在本文中，我们将尝试用MATLAB实现Romberg积分法，探索其优势和应用。
二、Romberg积分法原理Romberg积分法的基本原理是通过对梯形法则和辛普森法则进行逐步的细分和修正，以获得更加精确的数值积分结果。
假设我们需要求解函数 f(x) 在区间 [a, b] 上的定积分，那么Romberg积分法的步骤可以概括为以下几点：1. 将区间 [a, b] 均匀分成若干个子区间；2. 计算每个子区间上的梯形规则和辛普森规则的数值积分；3. 利用已知结果进行Richardson外推，修正数值积分的误差；4. 逐步增加子区间的细分，直到达到所需的精度要求。
三、MATLAB实现Romberg积分法我们可以使用MATLAB编程语言来实现Romberg积分法，以下是一个示例代码：function [I, R] = romberg(f, a, b, n)h = (b - a) ./ (2 .^ (0:n-1));R = zeros(n, n);R(1, 1) = (b - a) * (feval(f, a) + feval(f, b)) / 2;for j = 2:nsubtotal = 0;for i = 1:2^(j-2)subtotal = subtotal + feval(f, a + (2*i - 1) * h(j));endR(j, 1) = R(j-1, 1)/2 + h(j) * subtotal;endfor k = 2:nfor j = k:nR(j, k) = (4^(k-1) * R(j, k-1) - R(j-1, k-1)) / (4^(k-1) - 1); endendI = R(n, n);通过以上的MATLAB代码，我们可以轻松地实现Romberg积分法，对给定的函数和区间进行数值积分，并得到精确的积分结果。
MATLAB数值积分及算例
6.2.3 被积函数由一个表格定义
(要求积分,但是函数没有直接给出,只是自己在 做实验时得到的一组相关联的数据)
在MATLAB中，对由表格形式定义的函数关系的求定积分问 题用trapz(X,Y)函数。其中向量X,Y定义函数关系Y=f(X)。
例4 用trapz函数计算定积分。
命令如下：
X=1:0.01:2.5; Y=exp(-X); trapz(X,Y)
例2 求定积分:
x sin x
dx
0 (1 cos x cos x)
(1) 被积函数文件fx.m。
function f=fx(x) f=x.*sin(x)./(1+cos(x).*cos(x));
(2) 调用函数quad8求定积分。
I=quad8('fx',0,pi)
例3
分别用quad函数和quad8函数求定积分
global ki;ki=0; I=dblquad('fxy',-2,2,-1,1) ki
6.2 数值积分的实现方法
6.2.1 变步长辛普生法
基于变步长辛普生法，MATLAB给出了quad函数来 求定积分。该函数的调用格式为：
[I, n] = quad('fname', a, b, tol, trace) 其中fname是被积函数名。a和b分别是定积分的下 限 和 上 限 。 tol 用 来 控 制 积 分 精 度 ， 缺 省 时 取 tol=0.001。trace控制是否展现积分过程，若取非0则 展现积分过程，取0则不展现，缺省时取trace=0。 返回参数I即定积分值，n为被积函数的调用次数。
2.5 exdx
1
的近似值，并在相同的积分精度下，比较函数的调
Matlab数值积分程序集合
Matlab数值积分程序集合[图书馆+网络收集]近来学习数值积分,手头积累了不少程序,也拿来和各位朋友分享一下。
主要是来自数值积分教材和网络，基本的原理也就不打算多说了，随便搜索一下就可以得到，那就开始上代码了，呵呵，非原创，但是全部验证过，有疑问可以给我e-mail：1 梯形数值积分的MATLAB主程序function T=rctrap(fun,a,b,m)%fun 函数，a 积分上限 b积分下限 m 递归次数n=1;h=b-a; T=zeros(1,m+1); x=a;T(1)=h*(feval(fun,a)+feval(fun,b))/2;for i=1:mh=h/2; n=2*n; s=0;for k=1:n/2x=a+h*(2*k-1); s=s+feval(fun,x);endT(i+1)=T(i)/2+h*s;endT=T(1:m);e.g运行后屏幕显示精确值F s，用rctrap计算的递归值T和T与精确值F s的绝对误差w T>> ) exp((-x^.2./2)./(sqrt(2*pi)))T=rctrap(fun,0,pi/2,14), syms tfi=int(exp((-t^2)/2)/(sqrt(2*pi)),t,0, pi/2);Fs= double(fi), wT= double(abs(fi-T))fun =@(x)exp((-x^.2./2)./(sqrt(2*pi)))T =Columns 1 through 71.4168 1.3578 1.3313 1.3195 1.3142 1.3119 1.3109Columns 8 through 141.3105 1.3103 1.3102 1.3102 1.3101 1.3101 1.3101Fs =0.4419wT =Columns 1 through 70.9749 0.9159 0.8894 0.8776 0.8723 0.8700 0.8690Columns 8 through 140.8686 0.8684 0.8683 0.8683 0.8683 0.8682 0.8682>>2 复合辛普森（Simpson）数值积分的MATLAB主程序function y=comsimpson(fun,a,b,n)% fun 函数 a 积分上限 b积分下限 n 分割小区间数z1=feval (fun,a)+ feval (fun,b);m=n/2;h=(b-a)/(2*m); x=a;z2=0; z3=0; x2=0; x3=0;for k=2:2:2*mx2=x+k*h; z2= z2+2*feval (fun,x2);for k=3:2:2*mx3=x+k*h; z3= z3+4*feval (fun,x3);endy=(z1+z2+z3)*h/3;由于Matlab自带了 quad就是这个算法所以比较少自己编3 龙贝格数值积分的MATLAB主程序function [RT,R,wugu,h]=romberg(fun,a,b, wucha,m)%fun被积函数 a，b积分上下限 wucha两次相邻迭代绝对差值 m 龙贝格积分表最大行数%RT 龙贝格积分表 R 数值积分结果 wucha 误差估计 h 最小步长n=1;h=b-a; wugu=1; x=a;k=0; RT=zeros(4,4);RT(1,1)=h*(feval(fun,a)+feval(fun,b))/2;while((wugu>wucha)&(k<m)|(k<4))k=k+1; h=h/2; s=0;for j=1:nx=a+h*(2*j-1); s=s+feval(fun,x);endRT(k+1,1)= RT(k,1)/2+h*s; n=2*n;for i=1:kRT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1);endwugu=abs(RT(k+1,k)-RT(k+1,k+1));endR=RT(k+1,k+1);>> F=inline('1./(1+x)'); [RT,R,wugu,h]=romberg(F,0,1.5,1.e-8,13) syms xfi=int(1/(1+x),x,0,1.5); Fs=double(fi),wR=double(abs(fi-R)), wR1= wR - wuguRT =1.0500 0 0 00 00.9536 0.9214 0 00 00.9260 0.9168 0.9165 00 00.9187 0.9163 0.9163 0.91630 00.9169 0.9163 0.9163 0.9163 0.91630.9164 0.9163 0.9163 0.9163 0.9163 0.9163R =0.9163wugu =2.9436e-011h =0.0469Fs =0.9163wR =9.8007e-011wR1 =6.8571e-011>>6 复合梯形法function [I,step] = CombineTraprl(f,a,b,eps)%f 被积函数%a,b 积分上下限%eps 精度%I 积分结果%step 积分的子区间数if(nargin ==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f ),findsym(sym(f)),x1));endendI=I2;step=n;7 辛普森法function [I,step] = IntSimpson(f,a,b,type,eps)%type = 1 辛普森公式%type = 2 辛普森3/8公式%type = 3 复合辛普森公式if(type==3 && nargin==4)eps=1.0e-4; %缺省精度为0.0001endI=0;switch typecase 1,I=((b-a)/6)*(subs(sym(f),findsym(sym(f)),a)+...4*subs(sym(f),findsym(sym(f)),(a+b)/2)+...subs(sym(f),findsym(sym(f)),b));step=1;case 2,I=((b-a)/8)*(subs(sym(f),findsym(sym(f)),a)+...3*subs(sym(f),findsym(sym(f)),(2*a+b)/3)+ ...3*subs(sym(f),findsym(sym(f)),(a+2*b)/3)+subs(sym( f),findsym(sym(f)),b));step=1;case 3,n=2;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym (sym(f)),b))/h;while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/6)*(subs(sym(f),findsym(sym(f )),x)+...4*subs(sym(f),findsym(sym(f)), (x+x1)/2)+...subs(sym(f),findsym(sym(f)),x1 ));endendI=I2;step=n;end8 牛顿－科茨法function I = NewtonCotes(f,a,b,type)%type = 1 科茨公式%type = 2 牛顿－科茨六点公式%type = 3 牛顿－科茨七点公式I=0;switch typecase 1,I=((b-a)/90)*(7*subs(sym(f),findsym(sym(f)),a)+...32*subs(sym(f),findsym(sym(f)),(3*a+b)/4)+...12*subs(sym(f),findsym(sym(f)),(a+b)/2)+...32*subs(sym(f),findsym(sym(f)),(a+3*b)/4)+7*su bs(sym(f),findsym(sym(f)),b));case 2,I=((b-a)/288)*(19*subs(sym(f),findsym(sym(f)),a)+...75*subs(sym(f),findsym(sym(f)),(4*a+b)/5)+...50*subs(sym(f),findsym(sym(f)),(3*a+2*b)/5)+...50*subs(sym(f),findsym(sym(f)),(2*a+3*b)/5)+...75*subs(sym(f),findsym(sym(f)),(a+4*b)/5)+19*s ubs(sym(f),findsym(sym(f)),b));case 3,I=((b-a)/840)*(41*subs(sym(f),findsym(sym(f)),a)+...216*subs(sym(f),findsym(sym(f)),(5*a+b)/6)+...27*subs(sym(f),findsym(sym(f)),(2*a+b)/3)+...272*subs(sym(f),findsym(sym(f)),(a+b)/2)+...27*subs(sym(f),findsym(sym(f)),(a+2*b)/3)+...216*subs(sym(f),findsym(sym(f)),(a+5*b)/6)+41* subs(sym(f),findsym(sym(f)),b));end9 高斯公式function I = IntGauss(f,a,b,n,AK,XK)if(n<5 && nargin == 4)AK = 0;XK = 0;elseXK1=((b-a)/2)*XK+((a+b)/2);I=((b-a)/2)*sum(AK.*subs(sym(f),findsym(f),XK1));endta = (b-a)/2;tb = (a+b)/2;switch ncase 0,I=2*ta*subs(sym(f),findsym(sym(f)),tb);case 1,I=ta*(subs(sym(f),findsym(sym(f)),ta*0.5773503+tb)+...subs(sym(f),findsym(sym(f)),-ta*0.5773503+tb));case 2,I=ta*(0.55555556*subs(sym(f),findsym(sym(f)),ta*0.7745 967+tb)+...0.55555556*subs(sym(f),findsym(sym(f)),-ta*0.7 745967+tb)+...0.88888889*subs(sym(f),findsym(sym(f)),tb));case 3,I=ta*(0.3478548*subs(sym(f),findsym(sym(f)),ta*0.86113 63+tb)+...0.3478548*subs(sym(f),findsym(sym(f)),-ta*0.86 11363+tb)+...0.6521452*subs(sym(f),findsym(sym(f)),ta*0.339 8810+tb)...+0.6521452*subs(sym(f),findsym(sym(f)),-ta*0.3 398810+tb));case 4,I=ta*(0.2369269*subs(sym(f),findsym(sym(f)),ta*0.90617 93+tb)+...0.2369269*subs(sym(f),findsym(sym(f)),-ta*0.90 61793+tb)+...0.4786287*subs(sym(f),findsym(sym(f)),ta*0.538 4693+tb)...+0.4786287*subs(sym(f),findsym(sym(f)),-ta*0.5 384693+tb)+...0.5688889*subs(sym(f),findsym(sym(f)),tb)); end10 区间逐次分半梯形法function q=DblTraprl(f,a,A,b,B,m,n)if(m==1 && n==1) %梯形公式q=((B-b)*(A-a)/4)*(subs(sym(f),findsym(sym(f)),{a,b})+...subs(sym(f),findsym(sym(f)),{a,B})+...subs(sym(f),findsym(sym(f)),{A,b})+...subs(sym(f),findsym(sym(f)),{A,B}));else %复合梯形公式 C=4*ones(n+1,m+1);C(1,:)=2;C(:,1)=2;C(n+1,:)=2;C(:,m+1)=2;C(1,1)=1;C(1,m+1)=1;C(n+1,1)=1;C(n+1,m+1)=1; %C矩阵endF=zeros(n+1,m+1);q=0;for i=0:nfor j=0:mx=a+i*(A-a)/n;y=b+j*(B-b)/m;F(i+1,j+1)=subs(sym(f),findsym(sym(f)),{x,y});q=q+F(i+1,j+1)*C(i+1,j+1);endendq=((B-b)*(A-a)/4/m/n)*q;11 区间逐次分半布尔法function [I,step] = DDBuer(f,a,b,eps)if(nargin==3)eps=1.0e-4;end;n=1;h=b-a;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h/ 2;tol=1;while tol>epsn=n+1;h=h/4;I1=I2;I2=0;for i=0:(4^(n-1)-1)x=a+h*4*i;x1=x+h;x2=x1+h;x3=x2+h;x4=x3+h;I2=I2+(2*h/45)*(7*subs(sym(f),findsym(sym(f)),x)+...32*subs(sym(f),findsym(sym(f)),x1)+...12*subs(sym(f),findsym(sym(f)),x2)+...32*subs(sym(f),findsym(sym(f)),x3)+...7*subs(sym(f),findsym(sym(f)),x4));endtol=abs(I2-I1);endI=I2;step=n;10 自适应求积法function I=SmartSimpson(f,a,b,eps)if(nargin==3)eps=1.0e-4;end;e=5*eps;I=SubSmartSimpson(f,a,b,e);function q=SubSmartSimpson(f,a,b,eps)QA=IntSimpson(f,a,b,1,eps);QLeft=IntSimpson(f,a,(a+b)/2,1,eps);QRight=IntSimpson(f,(a+b)/2,b,1,eps);if(abs(QLeft+QRight-QA)<=eps)q=QA;elseq=SubSmartSimpson(f,a,(a+b)/2,eps)+SubSmartSimpson(f,(a+b)/2,b ,eps);ende.g>> SmartSimpson('x*sin(x)',0,1)ans =0.3011>>。
(完整版)数值积分及matlab实现
建立数值积分公式的途径比较多, 其中最常用的
有两种：
（1）由积分中值定理可知，对于连续函数f(x)，在
积分区间[a,b]内存在一点ξ，使得
分,因此将 选(x取) 为插值多项式, 这样f(x)的积分就
可以用其插值多项式的积分来近似代替
2.2 插值求积公式
设已知f(x)在节点 xk (k 0,1, , n) 有函数值 f (xk ) ,作n次拉格朗日插值多项式
式中
n
P(x) f (xk )lk (x)
k 0
lk (x)
n j0
b
n
f (x)dx
a
Ak f (xk )
k 0
为插值型求积公式的充要条件是公式
(
x)dx
时，则称求积公式为插值
设插值求积公式的余项为 R( f ) ,由插值余项定理得
R( f ) b f (x) P(x)dx b f (n1) ( ) (x)dx
a
a (n 1)!
其中 a, b
当f(x)是次数不高于n的多项式时，有 f (n1) (x) 0 R( f ) =0,求积公式(4)能成为准确的等式。由于闭区 间[a,b]上的连续函数可用多项式逼近，所以一个 求积公式能对多大次数的多项式f(x)成为准确等式， 是衡量该公式的精确程度的重要指标，为此给出以 下定义。
数值积分与微分
2009.4.22
数值积分和数值微分
1 引言 我们知道,若函数f(x)在区间[a,b]上连续且其原
函数为F(x),则可用Newton-Leibnitz公式
数值分析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。
matlab 数组积分
matlab 数组积分在MATLAB中，数值积分是常见的数值计算任务之一。
数值积分是对函数在给定区间上的积分值进行数值计算的过程。
在MATLAB中，有几种不同的方法可以用来进行数值积分。
一、MATLAB中的积分函数MATLAB提供了一些内置的函数，可以用来进行数值积分计算。
其中最常用的函数是`integral`函数。
`integral`函数可以用于一维和多维积分，可以使用固定步长或自适应步长算法。
下面是一个使用`integral`函数计算一维积分的示例：```matlabf = @(x) exp(-x^2); %定义需要积分的函数a = -1; %积分下限b = 1; %积分上限result = integral(f, a, b); %计算积分disp(result); %输出结果```在这个示例中，我们首先定义了需要积分的函数`f`，然后定义了积分的下限`a`和上限`b`。
然后我们使用`integral`函数来计算积分的值，并将结果存储在`result`变量中。
最后，我们使用`disp`函数来输出积分的结果。
除了`integral`函数，MATLAB还提供了其他一些积分函数，如`quad`、`quadl`、`quadgk`等。
这些函数提供了不同的积分算法和参数设置，可以根据具体的需求选择合适的函数进行数值积分计算。
二、积分方法在进行数值积分时，常用的方法包括：1.矩形法：将积分区间划分为若干个子区间，然后在每个子区间上选取某个点的函数值作为近似值。
这种方法简单易懂，但精度较低。
2.梯形法：将积分区间划分为若干个子区间，然后在每个子区间上通过线性插值得到函数的近似值，再对近似值进行积分。
这种方法比矩形法精度更高，但仍然有误差。
3.辛普森法：将积分区间划分为若干个子区间，然后在每个子区间上使用二次插值得到函数的近似值，再对近似值进行积分。
这种方法的精度比梯形法更高，但计算量也更大。
三、示例下面我们通过一个具体的示例来演示如何在MATLAB中进行数值积分计算。
如何在Matlab中进行数值积分和数值解
如何在Matlab中进行数值积分和数值解在数学和工程领域，数值积分和数值解是常见的技术手段，可以帮助我们求解复杂的数学问题和实际工程中的模型。
本文将介绍如何使用Matlab进行数值积分和数值解，以及一些注意事项和常用的方法。
一、数值积分数值积分是计算定积分的近似值的方法，可以通过数值逼近或数值插值来实现。
在Matlab中，有几种常用的函数可以用于数值积分，比如trapz、quad等。
1. trapz函数trapz函数是用梯形法则计算积分的函数。
它的使用方法是将要积分的函数作为输入的第一个参数，x轴上的点作为输入的第二个参数。
例如，要计算函数f(x)在区间[a, b]上的积分，可以使用以下代码：result = trapz(x, f(x));2. quad函数quad函数是使用自适应数值积分算法计算积分的函数。
它的使用方法是将要积分的函数作为输入的第一个参数，积分区间的下限和上限作为输入的第二个和第三个参数。
例如，要计算函数f(x)在区间[a, b]上的积分，可以使用以下代码：result = quad(@(x) f(x), a, b);二、数值解数值解是使用数值方法求解复杂的数学问题或实际工程中的模型的近似解。
在Matlab中，有几种常用的函数可以用于数值解，比如fsolve、ode45等。
1. fsolve函数fsolve函数是用于求解非线性方程组的函数。
它的使用方法是将非线性方程组表示为一个函数，然后将该函数作为输入的第一个参数。
例如，要求解方程组f(x) = 0，可以使用以下代码：x = fsolve(@(x) f(x), x0);其中x0是方程的初始猜测值。
2. ode45函数ode45函数是求解常微分方程初值问题的函数。
它的使用方法是将微分方程表示为一个函数，然后将该函数作为输入的第一个参数。
例如，要求解常微分方程dy/dx = f(x, y)，可以使用以下代码：[t, y] = ode45(@(t, y) f(t, y), tspan, y0);其中tspan是时间区间，y0是初始条件。
matlab中的微分方程的数值积分
MATLAB是一种流行的数学软件，用于解决各种数学问题，包括微分方程的数值积分。
微分方程是许多科学和工程问题的数学描述方式，通过数值积分可以得到微分方程的数值解。
本文将介绍在MATLAB中如何进行微分方程的数值积分，以及一些相关的技巧和注意事项。
一、MATLAB中微分方程的数值积分的基本方法1. 常微分方程的数值积分在MATLAB中，常微分方程的数值积分可以使用ode45函数来实现。
ode45是一种常用的数值积分函数，它使用4阶和5阶Runge-Kutta 方法来求解常微分方程。
用户只需要将微分方程表示为函数的形式，并且提供初值条件，ode45就可以自动进行数值积分，并得到微分方程的数值解。
2. 偏微分方程的数值积分对于偏微分方程的数值积分，在MATLAB中可以使用pdepe函数来实现。
pdepe可以求解具有定解条件的一维和二维偏微分方程，用户只需要提供偏微分方程的形式和边界条件，pdepe就可以进行数值积分，并得到偏微分方程的数值解。
二、在MATLAB中进行微分方程数值积分的注意事项1. 数值积分的精度和稳定性在进行微分方程的数值积分时，需要注意数值积分的精度和稳定性。
如果数值积分的精度不够，可能会导致数值解的误差过大；如果数值积分的稳定性差，可能会导致数值解发散。
在选择数值积分方法时，需要根据具体的微分方程来选择合适的数值积分方法，以保证数值解的精度和稳定性。
2. 初值条件的选择初值条件对微分方程的数值解有很大的影响，因此在进行微分方程的数值积分时，需要选择合适的初值条件。
通常可以通过对微分方程进行分析，或者通过试验求解来确定合适的初值条件。
3. 数值积分的时间步长在进行微分方程的数值积分时，需要选择合适的时间步长，以保证数值积分的稳定性和效率。
选择时间步长时，可以通过试验求解来确定合适的时间步长，以得到最优的数值解。
三、MATLAB中微分方程数值积分的实例以下通过一个简单的例子来演示在MATLAB中如何进行微分方程的数值积分。
数值积分matlab
数值积分matlab数值积分是数学中的一个重要分支，它是通过数值方法来计算函数的积分值。
在实际应用中，很多函数的积分无法通过解析方法求得，这时候就需要使用数值积分的方法来进行计算。
在本文中，我们将介绍如何使用Matlab进行数值积分。
Matlab是一种强大的数学计算软件，它提供了许多数值积分的函数，包括quad、quadl、quadgk等。
这些函数可以用来计算一维和多维函数的积分值。
我们来看一维函数的积分计算。
假设我们要计算函数f(x)在区间[a,b]上的积分值，可以使用quad函数来进行计算。
quad函数的语法格式为：[q,err] = quad(fun,a,b,tol)其中，fun是要计算的函数句柄，a和b是积分区间的上下限，tol 是误差容限。
函数quad会返回积分值q和误差err。
例如，我们要计算函数f(x)=x^2在区间[0,1]上的积分值，可以使用以下代码：fun = @(x) x.^2;[q,err] = quad(fun,0,1);运行以上代码，可以得到积分值q=0.3333和误差err=1.1102e-16。
接下来，我们来看多维函数的积分计算。
假设我们要计算函数f(x,y)在矩形区域[a,b]×[c,d]上的积分值，可以使用quad2d函数来进行计算。
quad2d函数的语法格式为：[q,err] = quad2d(fun,a,b,c,d)其中，fun是要计算的函数句柄，a、b、c和d是积分区间的上下限。
函数quad2d会返回积分值q和误差err。
例如，我们要计算函数f(x,y)=x^2+y^2在矩形区域[0,1]×[0,1]上的积分值，可以使用以下代码：fun = @(x,y) x.^2+y.^2;[q,err] = quad2d(fun,0,1,0,1);运行以上代码，可以得到积分值q=0.6667和误差err=1.1102e-16。
除了quad和quad2d函数外，Matlab还提供了许多其他的数值积分函数，如quadl、quadgk、integral等。
matlab数值积分实例
数值积分 1.求积分dx e x ⎰--112，在积分区间中，点与点之间的间隔取为0.1.解：（一）用MATBLE 编写复化梯形求积函数：function I=T_quad(x,y)n=length(x);m=length(y);if n ~=merrorendh=(x(n)-x(1))/(n-1);a=[1 2*ones(1,n-2) 1];I =h/2*sum(a.*y);输入：x=-1:0.1:1;y=exp(-x.^2);I =T_quad(x,y)运行得到：I =1.4924(二)用MATBL 编写复化Simpson 求积函数：function I=S_quad(x,y);n=length(x);m=length(y);if n ~=merrorendif rem(n-1,2) ~=0I=T_quad(x,y);return;endN=(n-1)/2;h=(x(n)-x(1))/N;a=zeros(1,n);for k=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4;a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);输入：x= -1:0.1:1;y=exp(-x.^2);I= S_quad(x,y)运行得到：I =1.4936（三）用MATBL 编写复化Cotes 求积函数：function I=C_quad(x,y);n=length(x);m=length(y);if n ~=merrorendif rem(n-1,4) ~=0I=S_quad(x,y);returnendN=(n-1)/4;h=(x(n)-x(1))/N;a=zeros(1,n);for k=1:Na(4*k-3)=a(4*k-3)+7;a(4*k-2)=a(4*k-2)+32;a(4*k-1)=a(4*k-1)+12;a(4*k)=a(4*k)+32;a(4*k+1)=a(4*k+1)+7;endI=h/90*sum(a.*y);输入：x= -1:0.1:1;y=exp(-x.^2);I= C_quad(x,y)运行得：I =1.4936（四）利用trapz()函数，采用复化梯形公式求积分输入：x=-1:0.1:1;y=exp(-x.^2);I=trapz(x,y)输出：I =1.49242. 求积分dx e x ⎰--112，取精度要求510-=ε （一）用MATLAB 编写自适应步长的梯形公式function I=T_quad_iter(fun,a,b,ep)if nargin<4 ep=1e-5;endN=1;h=b-aT=h/2*(feval(fun,a)+feval(fun,b));while 1h=h/2;I=T/2;for k=1:NI=I+h*feval(fun,a+(2*k-1)*h);endif abs(I-T)<ep breakendN=2*N;T=IEnd输入：I=T_quad_iter(@(x)exp(-x.^2),-1,1)输出：I =1.4936（二）用MATLAB编写自适应步长的Simpson公式function I=S_quad_iter(fun,a,b,ep)if nargin<4 ep=1e-5endN=1;h=b-aT1=h/2*(feval(fun,a)+feval(fun,b));S0=T1;while 1h=h/2;T2=T1/2;for k=1:NT2=T2+h*feval(fun,a+(2*k-1)*h)endI=(4*T2-T1)/3;if abs(I-S0)<ep breakendN=2*N;T1=T2;S0=IEnd输入：I=S_quad_iter(@(x) exp(-x.^2),-1,-1)输出：I =1.4936（三）用quad（）函数，采用自适应步长的Simposon求积分输入：I=quad(@(x) exp(-x.^2) ,-1,1)输出：I =1.4936。
应用MATLAB进行实验数据的数值积分
实验数据的数值积分问题用不同温度下固体铅的定压热容数据，求298K 时的绝对熵。
数据见下表：由热力学方法可知，绝对熵与定压热容有如下关系：dT TC S TPT ⎰=0采用MATLAB 来计算，这里将最小二乘B 样条拟合法和三次样条插值函数计算作比较，命令行程序如下：function Entropy_Int clear all clc% 读入数据T = [5. 10. 15. 20. 25. 30. 50. 70. 100. 150. 200. 250. 298.];Cp = [0.305 2.8 7.0 10.8 14.1 16.5 21.4 23.3 24.5 25.4 25.8 26.2 26.5]; CpT = Cp./T;T = [0 T]; CpT = [0 CpT];% 用最小二乘B 样条拟合法和三次样条插值函数计算 knots = 3; K = 3; % 三次B 样条 sp = spapi(knots,K,T,CpT);cs = csapi(T,CpT); % 生成三次样条插值函数cs % 绘制浓度拟合曲线Ti = linspace(T(1),T(end),299); CpT_B = fnval(sp,Ti); CpT_c = fnval(cs,Ti);plot(T,CpT,'ro',Ti,CpT_B,'b-',Ti,CpT_c,'Bl-.') xlabel('T') ylabel('Cp/T')legend('实验值','B 样条拟合','立方样条拟合') % 进行数值积分 pp = fnint(cs); s = fnval(pp,T)TC p /Ts =0 0.0433 0.8733 2.7820 5.3352 8.1104 10.9071 20.7236 28.2599 36.8085 46.9207 54.2919 60.0880 64.7233在解决本问题时，需首先要考虑积分下限，即0=T 时的情形，因此时的C P 为0/0。
[转载]matlab数值积分实现
![[转载]matlab数值积分实现](https://img.taocdn.com/s3/m/01783b3d5b8102d276a20029bd64783e09127d7d.png)
[转载]matlab数值积分实现原⽂地址：作者：最近做有关加速度的数据处理，需要把加速度积分成位移，⽹上找了找相关资料，发现做这个并不多，把最近做的总结⼀下吧！积分操作主要有两种⽅法：时域积分和频域积分，积分中常见的问题就是会产⽣⼆次趋势。
关于积分的⽅法，在国外⼀个论坛上有⼈提出了如下说法，供参考。
Double integration of raw acceleration data is a pretty poor estimate for displacement. The reason is that at each integration, you are compounding the noise in the data.If you are dead set on working in the time-domain, the best results come from the following steps.1. Remove the mean from your sample (now have zero-mean sample)2. Integrate once to get velocity using some rule (trapezoidal, etc.)3. Remove the mean from the velocity4. Integrate again to get displacement.5. Remove the mean. Note, if you plot this, you will see drift over time.6. To eliminate (some to most) of the drift (trend), use a least squares fit (high degree depending on data) to determine polynomial coefficients.7. Remove the least squares polynomial function from your data.A much better way to get displacement from acceleration data is to work in the frequency domain. To do this, follow these steps...1. Remove the mean from the accel. data2. Take the Fourier transform (FFT) of the accel. data.3. Convert the transformed accel. data to displacement data by dividing each element by -omega^2, where omega is the frequency band.4. Now take the inverse FFT to get back to the time-domain and scale your result.This will give you a much better estimate of displacement.说到底就是频域积分要⽐时域积分效果更好，实际测试也发现如此。
matlab中求积分的命令
matlab中求积分的命令求积分是数学中的一个重要概念，也是数学分析中的基础内容。
在MATLAB中，我们可以使用一些特定的命令来实现对函数的积分计算，从而得到函数的解析式或数值结果。
本文将介绍一些常用的MATLAB求积分命令，并探讨其在实际问题中的应用。
一、MATLAB中的求积分命令在MATLAB中，求积分的命令主要有两种：符号积分和数值积分。
下面分别介绍这两种求积分的命令及其使用方法。
1. 符号积分命令符号积分是指对给定的函数进行解析求积分，得到一个含有未知常数的解析式。
在MATLAB中，可以使用符号积分命令'int'来进行符号积分的计算。
其基本语法为：int(f, x) 或 int(f, x, a, b)其中，f表示被积函数，x表示积分变量，a和b表示积分区间的上下限。
例如，要对函数f(x) = x^2进行符号积分，可以使用以下命令：syms xf = x^2;F = int(f, x)这样，MATLAB将输出函数F(x) = (1/3)x^3，即f(x)的积分结果。
2. 数值积分命令数值积分是指对给定的函数进行数值近似求积分，得到一个数值结果。
在MATLAB中，可以使用数值积分命令'integral'来进行数值积分的计算。
其基本语法为：Q = integral(fun, a, b)其中，fun表示被积函数的函数句柄，a和b表示积分区间的上下限。
例如，要对函数f(x) = exp(-x^2)进行数值积分，可以使用以下命令：f = @(x) exp(-x^2);Q = integral(f, -inf, inf)这样，MATLAB将输出数值结果Q，即f(x)的积分值。
二、MATLAB求积分命令的应用MATLAB中的求积分命令在工程和科学计算中有着广泛的应用。
下面将介绍两个实际问题的求解过程，以展示这些命令的应用。
1. 求解概率密度函数的积分概率密度函数是统计学中的一个重要概念，用于描述随机变量的概率分布。
三个参数matlab程序,用matlab求定积分的三个实例代码
三个参数matlab程序,⽤matlab求定积分的三个实例代码⼀、符号积分符号积分由函数int来实现。
该函数的⼀般调⽤格式为：int(s)：没有指定积分变量和积分阶数时，系统按findsym函数指⽰的默认变量对被积函数或符号表达式s求不定积分；int(s,v)：以v为⾃变量，对被积函数或符号表达式s求不定积分；int(s,v,a,b)：求定积分运算。
a,b分别表⽰定积分的下限和上限。
该函数求被积函数在区间[a,b]上的定积分。
a和b可以是两个具体的数，也可以是⼀个符号表达式，还可以是⽆穷(inf)。
当函数f关于变量x在闭区间[a,b]上可积时，函数返回⼀个定积分结果。
当a,b中有⼀个是inf 时，函数返回⼀个⼴义积分。
当a,b中有⼀个符号表达式时，函数返回⼀个符号函数。
例：求函数x^2+y^2+z^2的三重积分。
内积分上下限都是函数，对z积分下限是sqrt(x*y)，积分上限是x^2*y；对y积分下限是sqrt(x)，积分上限是x^2；对x的积分下限1，上限是2，求解如下：>>syms x y z %定义符号变量>>F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) %注意定积分的书写格式F2 =1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4) %给出有理数解>>VF2=vpa(F2) %给出默认精度的数值解VF2 =224.92153573331143159790710032805⼆、数值积分1.数值积分基本原理求解定积分的数值⽅法多种多样，如简单的梯形法、⾟普⽣(Simpson)法、⽜顿－柯特斯(Newton-Cotes)法等都是经常采⽤的⽅法。
它们的基本思想都是将整个积分区间[a,b]分成n个⼦区间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。
matlab数值计算代码
matlab数值计算代码Matlab是一种强大的数值计算软件，广泛应用于科学研究、工程设计等领域。
在Matlab中，我们可以使用代码来进行各种数值计算，包括数值积分、数值求解方程、数值解微分方程等。
本文将介绍一些常见的数值计算代码，并说明其原理和应用。
一、数值积分数值积分是利用数值方法求解定积分的过程。
在Matlab中，我们可以使用simpson函数或trapz函数进行数值积分计算。
这两个函数分别采用辛普森公式和梯形公式进行数值积分近似。
例如，下面的代码使用simpson函数计算函数f(x)在区间[a,b]上的定积分：```matlaba = 0;b = 1;n = 100;x = linspace(a, b, n);y = f(x);integral = simpson(y, x);```其中，a和b分别是积分区间的上下限，n是划分区间的个数，x是划分后的区间点，y是函数在各个区间点处的函数值，integral是计算得到的定积分值。
二、数值求解方程数值求解方程是指利用数值方法求解方程的近似解。
在Matlab中，我们可以使用fzero函数或fsolve函数进行数值求解方程。
这两个函数采用不同的求解算法，可以用于求解单变量方程或多变量方程。
例如，下面的代码使用fzero函数求解方程f(x)=0：```matlabx0 = 0;x = fzero(@f, x0);```其中，x0是求解初始值，@f是函数句柄，表示要求解的方程。
x 是求解得到的近似解。
三、数值解微分方程数值解微分方程是指利用数值方法求解微分方程的近似解。
在Matlab中，我们可以使用ode45函数或ode23函数进行数值解微分方程。
这两个函数采用不同的数值方法，可以用于求解常微分方程或偏微分方程。
例如，下面的代码使用ode45函数求解常微分方程dy/dx=f(x,y)：```matlabx0 = 0;y0 = 1;xspan = [0, 1];[t, y] = ode45(@f, xspan, y0);```其中，x0和y0分别是初始条件，xspan是求解区间，@f是函数句柄，表示要求解的微分方程。
【MATLAB】实验五：数值微积分与方程数值求解
实验五 数值微积分与方程数值求解一、实验目的1. 掌握求数值导数和数值积分的方法。
2. 掌握代数方程数值求解的方法。
3. 掌握常微分方程数值求解的方法。
二、实验内容要求：命令手工 ( )输入1. 求函数在指定点的数值导数。
232()123,1,2,3026x x x f x x x x x==2. 用数值方法求定积分。
(1) 210I π=⎰的近似值。
(2) 2220ln(1)1x I dt xπ+=+⎰3. 分别用三种不同的数值方法解线性方程组。
6525494133422139211x y z u x y z u x y z u x y u +-+=-⎧⎪-+-=⎪⎨++-=⎪⎪-+=⎩4. 求非齐次线性方程组的通解。
1234123412342736352249472x x x x x x x x x x x x +++=⎧⎪+++=⎨⎪+++=⎩解：先建立M 函数文件,然后命令窗口中写命令。
121/119/112/115/111/1110/11100010X k k --⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=++⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦，其中12,k k 为任意常数。
5. 求代数方程的数值解。
(1) 3x +sin x -e x =0在x 0=1.5附近的根。
(2) 在给定的初值x 0=1，y 0=1，z 0=1下，求方程组的数值解。
23sin ln 70321050y x y z x z x y z ⎧++-=⎪+-+=⎨⎪++-=⎩ans =1289/6826. 求函数在指定区间的极值。
(1) 3cos log ()xx x x x f x e ++=在(0,1)内的最小值。
(2) 33212112122(,)2410f x x x x x x x x =+-+在[0,0]附近的最小值点和最小值。
（以下选作题，是微分方程的数值解）7. 求微分方程的数值解。
x 在[1.0e-9,20]2250(0)0'(0)0xd y dy y dx dx y y ⎧-+=⎪⎪⎪=⎨⎪=⎪⎪⎩解：M 文件：运行结果：8. 求微分方程组的数值解，并绘制解的曲线。
matlab 积分数值求解程序
一、简介Matlab是一种用于数学计算、数据分析和可视化的高级编程语言和环境。
其中，积分数值求解是Matlab中常用的功能之一，能够用于求解复杂的数学积分问题。
本文将介绍如何使用Matlab的积分数值求解程序，以及其基本原理和用法。
二、Matlab积分数值求解程序的基本原理Matlab中的积分数值求解程序主要基于数值积分的基本原理，即将一个函数在某个区间上的积分近似表示为若干个小区间上函数值的加权和。
常见的数值积分方法包括梯形法则、辛普森法则、龙贝格积分法等。
这些数值积分方法可以用于求解一维和多维函数的数值积分，并且在Matlab中都有相应的函数实现。
三、Matlab积分数值求解程序的使用方法1. 使用内置函数Matlab中提供了多个内置的积分数值求解函数，例如quad、quadl、quadgk等，这些函数可以用于对一维函数的数值积分。
用户只需输入被积函数、积分区间等参数即可得到积分的数值近似解。
2. 自定义函数除了使用内置函数，用户还可以通过自定义函数的方式来实现积分数值求解。
通过编写自己的积分数值求解算法，用户可以更灵活地对复杂的积分问题进行求解。
Matlab还提供了丰富的数学函数和工具，可以用于对被积函数进行求导、符号计算等，从而辅助完成数值积分的求解过程。
四、Matlab积分数值求解程序的应用实例1. 一维函数的数值积分假设有一个一维函数f(x)=x^2，在区间[0,1]上的数值积分。
可以使用Matlab的内置函数quad来求解这个积分，代码如下：```f = (x) x^2;q = quad(f, 0, 1);disp(['The value of the integral is ', num2str(q)]);```上述代码中，使用了匿名函数来定义被积函数f(x)=x^2，并将其作为参数传递给quad函数。
使用disp函数输出积分的数值近似解。
2. 多维函数的数值积分对于多维函数的数值积分，Matlab同样提供了相应的函数实现。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
（4）四阶亚当姆斯预估校正法：
clear all close all a=[0,1,0;0,0,1;-22.06,-27,-10]; b=[0;0;40.6]; u=1; for k=1:4 if (k==1) h=0.01; elseif (k==2) h=0.1; elseif (k==3) h=0.105; else (k==4) h=0.21; end x=[0;0;0]; y=x'; t(1)=0; for i=1:3 t(i+1)=h*i; k1=a*x+b*u; k2=a*(x+h/2*k1)+b*u; k3=a*(x+h/2*k2)+b*u; k4=a*(x+h*k3)+b*u; x=x+h/6*(k1+2*k2+2*k3+k4); y=[y;x']; y1=y'; end for i=4:1000 t(i+1)=h*i;
（3）给出仿真结果曲线图。
欧拉法：
改进欧拉法：
四阶龙格库塔法：
四阶亚当姆斯预估校正法：
实验结果得到四阶亚当姆斯预估校正法的临界步长约为 0.21。
5、实验结果及分析：四种方法的精度比较、稳定性比较，步长 稳定区是否符合实验前的理论分析值。其他需要总结的问题。
精度比较： 从保存的数据文件来看， 四阶龙格库塔法和四阶亚当姆斯预估校正法的精度 较高，大约数量级在 104 ，改进欧拉法次之，在 103 数量级上，欧拉法最低，在 102 数量 级上，与理论的结果相符程度较好。 稳定性比较：从实验结果看，四种方法在步长一样且较小的情况下，都能较快的达到稳 定。而且，同一种方法随着步长的增加，达到稳态的时间会增加。但是由于四种方法的临界 步长不相同，因而随着步长的增加，达到稳态的时间有所不同。 步长稳定区：从实验仿真的结果可以看出，实际的临界步长符合实验前理论的分析值。 而且四阶龙格库塔法的临界步长最大，四阶亚当姆斯预估校正法的临界步长最小。 本实验中四阶亚当姆斯预估校正法不能自启动， 需要用四阶龙格库塔法计算其起始的值。
（2）改进欧拉法：
clear all close all
a=[0,1,0;0,0,1;-22.06,-27,-10]; b=[0;0;40.6]; u=1; for k=1:4 if (k==1) h=0.01; elseif (k==2) h=0.1; elseif (k==3) h=0.16; else (k==4) h=0.321; end x=[0;0;0]; y=x'; t(1)=0; subplot(2,2,k) for i=1:1000 t(i+1)=h*i; x1=x+(a*x+b*u)*h; x=x+[(a*x+b*u)+(a*x1+b*u)]*h/2; y=[y;x']; end t1=0:0.01:10; y1=1.84-4.95*t1.*exp(-1.88*t1)-1.5*exp(-1.88*t1)-0.34*exp(-6.24*t 1); plot(t,y) hold on plot(t1,y1,'k-.','linewidth',2) xlabel('time(sec)') ylabel('·ùÖµy') legend('y','doty','ddoty','y-shuzhi') if (k==1) title('eulerp h=0.01') axis([0,10,-1.5,3.5]) elseif (k==2) title('eulerp h=0.1') axis([0,10,-1.5,3]) elseif (k==3) title('eulerp h=0.16 axis([0,10,-1.5,2]) else (k==4) title('eulerp h=0.321 临界') axis([0,50,-20,5]) 0.5*临界')
6、收获。
（1）这次试验我了解到连续系统仿真不同的数值积分法在精度上，步长稳定区以及稳 定性上有了较为直观的理解， 可以基本上掌握了不同数值积分法的选取在稳态精度上的不同。 （2）本次试验在 matlab 编程上有了较高的提高，在编程过程中，对 matlab 编程会出现 的问题以及解决方法有了基本上的掌握，收获较大。
；
（3）理论分析：计算系统特征值。结合系统特征值，对四种方法的稳定性进行分析，确定 每种方法步长的取值范围，即 h 临界。
1.8680 用 matlab 计算系统的特征值为： D 1.8928 ，可见系统的特征值有三个，我们取 6.2392
最大的特征值为： = 6.2392 。由四种方法的稳定区域可得： 方法 欧拉法 改进欧拉法 四阶龙格库塔法 四阶亚当姆斯法-显式 四阶亚当姆斯法-隐式 稳定区域 步长 h 临界步长 0.3206 0.3206 0.4456 0.048 0.2944
end end
（3）四阶龙哥库塔法：
clear all close all a=[0,1,0;0,0,1;-22.06,-27,-10]; b=[0;0;40.6]; u=1; for k=1:4 if (k==1) h=0.01; elseif (k==2) h=0.1; elseif (k==3) h=0.223; else (k==4) h=0.447; end x=[0;0;0]; y=x'; t(1)=0; subplot(2,2,k) for i=1:1000 t(i+1)=h*i; k1=a*x+b*u; k2=a*(x+h/2*k1)+b*u; k3=a*(x+h/2*k2)+b*u; k4=a*(x+h*k3)+b*u; x=x+h/6*(k1+2*k2+2*k3+k4); y=[y;x']; end t1=0:0.01:10; y1=1.84-4.95*t1.*exp(-1.88*t1)-1.5*exp(-1.88*t1)-0.34*exp(-6.24*t 1); plot(t,y) hold on plot(t1,y1,'k-.','linewidth',2) xlabel('time(sec)') ylabel('·ùÖµy') legend('y','doty', if (k==1) title('rk h=0.01') 'ddoty','y-shuzhi')
P yn yn f (tn , yn )h 1 改进欧拉法： ； h C P P yn 1 yn [ f (tn , yn ) f (tn 1 , yn 1 )] 2
k1 f (tn , yn ) k2 f (tn h , yn h k1 ) 2 2 h h 四阶龙哥库塔法： k3 f (tn , yn k2 ) ； 2 2 k4 f (tn h, yn hk3 ) yn 1 yn h (k1 2k2 2k3 k4 ) 6
(2,0)
(2,0) (2.78,0)
0 h 0.3206
0 h 0.3206 0 h 0.4456
0 h 0.048 0 h 0.2944
( 310 , 0)
( 90 49 , 0)
由于四阶亚当姆斯预估校正法综合受到显式与隐式的影响，其临界步长应在 0.048~0.2944 之间，具体值需用实验验证出来。
h P y y (55 f n 59 f n 1 37 f n 2 9 f n 3 ) n 1 n 24 四阶亚当姆斯预估校正法： y C y h 9 f P 19 f 5 f f n 1 n n 1 n n 1 n2 24
此处 U 为单位阶跃信号。
1 0 0 0 0 1 ，b 设： a 0 0 。 22.06 27 10 40.6
（2）分别写出四种方法的计算公式； 欧拉法： yn 1 yn f (tn , yn )h ；
控制系统仿真实验
姓名： 王体强 班 级 ： 120324 学 号 ： 12031190 指导教师： 王卫红 时 间 ： 2015 年 4 月 13 日
控制系统仿真实验
1、实目的：
进一步掌握数值积分法； 进一步掌握 MATLAB 软件的使用方法。
2、实验设备：
数字计算机，MATLAB 软件
3、实验预备：
（1）将传递函数化为一阶微分方程组（即状态方程） ； 设： y1 y ， y2 y ， y3 y
1 0 y1 0 y1 0 y 0 0 1 2 y2 0 U 22.06 27 10 40.6 y3 y3
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*临界')