数值积分法是求定积分近似值的数值方法。 也就是说，用被积函数的有限个采样值的离散或加权平均近似值代替定积分的值。
在求某个函数的定积分时，往往很难用初等函数表示被积函数的原始函数，而且很多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对于这类函数的定积分也无法用不定积分法求解。 在微积分学方面做出了杰出贡献的I.xhddp、l .欧拉、C.F .紧致发带、拉格朗日等数学大师，在数值积分这一领域做出了各自的贡献，为这一分支奠定了理论基础。
数值积分法也是计算机模拟中常用的方法。 已知函数微分方程时，求函数下一时刻的值，以欧拉法、梯形法、龙格库塔法为主。
欧拉法，这些方法中精度最低，程序比较简单。 欧拉法的公式可以写如下。
用欧拉法近似如下。
y'=f(t，y ) )。
这里，y’(k )函数y ) x )在k个时刻的导数，h是积分的步骤(也称为采样周期)。
欧拉法的主要想法是用当前时刻的值y(k )当前时刻y’) k ) h )积分步骤)近似代替y ) k1 )时刻的值。 这种方法的主要误差来源于y’(t )在h区间的变化，这里用y ) k )代替。
梯形法是基于欧拉法改进的算法，也称为改进的欧拉法。
用欧拉法求出y(k1 )后，求出y ) k1 )下的导数f ) k1，y ) k1 ) )。 用求梯形面积的方法求解该积分的大小，可以得到以上公式。
ymdsj(runge-kutta )法是工程上广泛应用的高精度单步算法，用于数值求解微分方程。 该算法精度高，采取措施抑制误差，因此其实现原理也很复杂。 这里主要介绍工程中常用的四阶ymdsj法。
该算法基于数学辅助。 一阶精度的拉格朗日中值定理包括：
微分方程时：
y'=f(x，y ) ) )。
y(I1 )=y (I ) h*K1
k1=f(Xi，yi ) ) ) ) ) ) ) ) ) )。
如果将点xi处斜率的近似值K1和右端点xi 1处的斜率K2的算术平均值设为平均斜率K*的近似值，则得到2次精度得到改善的拉格朗日中值定理：
y(I1 )=y (I ) [h* ) k1k2)/2]
k1=f(Xi，yi ) ) ) ) ) ) ) ) ) )。
(K2=F(x ) I ) h，y ) I ) h*K1 ) ) ) ) ) ) ) ) )
依次地，在区间[xi、xi 1]内，多估计几个点处梯度值K1、K2、……Km，将它们的加权平均作为平均梯度K*的近似值，就能够构筑精度明显高的高次计算式。 通过数学推导和求解，得到4阶龙格-库塔算法。 也就是说，在工程中应用广泛的经典龙格-库塔算法：
y(I1 )=y (I ) h* ) K12*K2*K3K4)/6
k1=f(x ) I )，y (I ) )
(K2=F(x ) I ) h/2，y ) I ) h*K1/2 ) ) ) ) ) ) )
(k3=f(x ) I ) h/2，y ) I ) h*K2/2 ) ) ) ) ) ) )
(k4=f(x ) I ) h，y ) I ) h*K3 ) ) ) ) ) ) ) ) )
k1是时间段开始时倾斜
k2是时间带中点的倾斜度，根据欧拉法采用倾斜度k1来决定y在点tn h/2处的值
k3也是中点的斜率，但这次采用斜率k2来决定y值；
k4是时间段终点的斜率，其y值由k3决定。
以上三种方法，ymdsj法精度最高，改进的梯形法次之，欧拉法精度最低。 关于精度，引入了截止误差的概念。
截止误差是指近似解与其真值之差。 将数值积分对的真值进行kdyd展开。 欧拉法在数值积分时，取kdyd级数的第一项，梯形法取kdyd级数的前两项，四阶ymdsj法取kdyd级数的前三项。 所以他们的截止误差是kdyd级数的其余部分。
数值积分
法在控制领域主要应用于对状态空间的输出进行求解。下面是我做过的一个matlab仿真已知开环传递函数
5 s + 100
-------------------------------------------
s^4 +8 s^3 + 31.99 s^2 + 75.21 s
给定单位阶跃，求解闭环函数的输出曲线。
首先我们要把传递函数转换成状态方程。在
主M文件：
clear
t=15;h=0.1;
num = [0 0 0 5 100];
den = [1 8 31.99 75.21 0];
den = den + num;
Tsys1 =tf(num,den);
[A,B,C,D] =tf2ss(num,den);
x0 = [0 0 0 0]';
u = 1;
i=0:h:t+h;
Y0 = EEluer(h,t,A,B,C,D,x0,u);
Y1 = SecondRunge(h,t,A,B,C,D,x0,u);
Y2 = FourRunge(h,t,A,B,C,D,x0,u);
plot(i,Y0,i,Y1,'-*',i,Y2,'g-');
legend('euler','SecondRunge','FourRunge');
四阶龙格库塔法M function：
function save=FourRunge(h,t,A,B,C,D,x0,u)
count = 1;
y0 =C * x0 + D * u;
len = size(y0,1);
wide=length(0:h:t);
save=zeros(wide,len);
save(count,len) = y0;
for i=0:h:t
k1 = h * My_function(A,B,x0,u) ;
k2 = h * My_function(A,B,x0+0.5*k1,u);
k3 = h * My_function(A,B,x0+0.5*k2,u);
k4 = h * My_function(A,B,x0+k3,u);
x1 = x0 + (1/6)*(k1 + 2*k2 +2*k3 +k4);
y0 = C * x1 + D * u;
count = count + 1;
save(count,len) = y0;
x0 = x1;
end
end
欧拉法 M function ：
function save=EEluer(h,t,A,B,C,D,x0,u)
count = 1;
y0 =C * x0 + D * u;
len = size(y0,1);
wide=length(0:h:t);
save=zeros(wide,len);
save(count,len) = y0;
for i=0:h:t
x1 = h * My_function(A,B,x0,u) +x0;
y0 = C * x1 + D * u;
count = count + 1;
save(count,len) = y0;
x0 = x1;
end
end
梯形法 M function：
function save=SecondRunge(h,t,A,B,C,D,x0,u)
count = 1;
y0 =C * x0 + D * u;
len = size(y0,1);
wide=length(0:h:t);
save=zeros(wide,len);
save(count,len) = y0;
for i=0:h:t
k1 = h * My_function(A,B,x0,u) ;
k2 = h * My_function(A,B,x0+k1,u);
x1 = x0 + (1/2)*(k1 + k2);
y0 = C * x1 + D * u;
count = count + 1;
save(count,len) = y0;
x0 = x1;
end
end
微分方程 M function
function [ yy ] = My_function(A,B,x0,u)
yy = A*x0 + B*u;
end
实验结果
图中振荡最大的是欧拉法，梯形法和四阶龙格库塔法精度差不多。四阶龙格库塔法精度最高。