几种数字仿真的物理意义与代码实现

任何的系统都包括输入与输出，小到一片最小的代码片，大道航天工业系统，都是由最基本的输入输出与中间环节构成。如何处理中间的环节就显得异常重要。现如今的系统大多为数字仿真系统，也就是大多使用计算机来处理输入与反馈信号。也信号又分为连续信号与离散信号，今天我们来谈一谈利用计算机仿真来处理信号的几种典型的方法及实现的MATLAB代码。

（1）数值积分法：

我们知道，从数学的角度上考虑，y(k+1)=y(k)+∫f(t,y(t))dt这样得来的确实是最优解，如果一直这样计算下去，得到的永远都是最精确地值。但是我们的计算机却无法做到这样的计算，因为两时刻之间的斜率为变化值。于是我们利用简化的方式去代替积分，只要采样时刻取得足够小，我们就可以理所当然的人为t（k）时刻与t（k+1）时刻之间是一条直线，于是我们就认为两者之间的斜率就是不变的值。这样看来就是好办多了只要曲积分为t(k)时刻的斜率，就解决了。下面为实例实现代码：

已知：y’ = -y; y(0) = 1;

function y = euler(y0,h,f)
y = y0 + h*f;
%代码如下 ：
h = 0.1;        %步长确定了精度；
y0 = 1;         %确定初值
t = 0:h:1;
n = length(t);
numy = zeros(1,n);
f = -y0;
numy(1) = 0;
for i = 2:n
numy(i) = euler(y0,h,f);
y0 = numy(i);
f = -y0;
end
realy = exp(-t);
plot(t,realy,t,numy,'ro-')
legend('真值'，'欧拉')
<注意：有的版本内euler函已经是新版本MATLAB的内部构造函数，无法再次定义使用，并且更改函数M文件是要保证函数名与被调用的函数名相同，不然代码会无法运行，你也无法调试出来错误！>

结果如下：


可见，这种方法仿真的结果随着时间的推移误差会越来越大，并且精度会随着步长的变大而差，但是这种方法提供了最基本的思想，所以一直是教学的典例。

（二）梯形法：y(k+1)=y(k)+(1/2)*h[f(tk,yk) + f(t(k+1),y(k+1))]，以上·公式可以看出一个问题，就是我要得到y（k+1）的值，利用梯形法，就必须先要的到y（k+1）的值，似乎变成了鸡生蛋的问题，这里我们利用欧拉方式先进行预测，然后结合梯形法，这样一来，梯形法也就是一种改进的欧拉方法，只是提高了精度。下面看看他是如何提高精度的，我们利用泰勒展开将t（k+1）时刻的值利用泰勒展开得到了如下的值：y（t+1） = y(k) + hy’(k) + (1/2!)h^2y”(t)·········

我们利用欧拉公式，得到了前两项，其舍去了2次方及其以后的部分，因此其截断误差也就是O（h^2）,相比下，梯形法虽然只是改变了其去斜率的方式，但是截断误差却提高了一个等级。

（三）龙哥-库塔法：

我们利用提梯形法提高了精度，但是如果仿真的时间过长，误差还是会一步一步的变差，提高精度的精髓在于如何设计积分中的斜率，于是引入龙哥-库塔法：

梯形法中我们去斜率我y（k）与y（k+1）的斜率的平均，也就是取得二者同样的权重，现在用同样的方法去两点之间的点，只要给点的斜率赋合理的权重，就可以很好的提高截断误差的阶次即精度。那么如何巧合的得提高截断误差？这里不得不佩服龙格的逆行思维，他采用倒推方式，先设置每个点的斜率的加权为一个未知数，带着未知数去推导，这样一来，将得到的y（k+1）时刻的值与泰勒展开进行比较，一一对比，从而得到了每一阶次所应该设定的权重。我们最多应用的也就是四阶龙格-库塔法。其取点的权重正好是（1 2 2 1），而二阶的龙格-库塔也正好与梯形法对应上。其实现代码：

function dy = ff(y0)
dy = -y0;
function funge4(y0,h)
k1 = h*ff(y0);
k2 = h*ff(y0 + k1/2);
k3 = h*ff(y0 + k3/2);
k4 = h*ff(y0 + k3);
yk_1 = y0 + (k1 + 2*k2 + 2*k3　+ k4)/6;
%下面为实现的main函数：
y0 = 1;
h = 0.1;
t = 0:h:1;
n = length(t);
runge4_y = zeros(1,n);
runge4_y(1) = 1;
for i = 2:n
runge4_y(i) = runge4(y0,h);
y0 = runge4_y(1);
end
realy = exp(-t);
plot(t,realy,t,runge4_y,'-o')

结果如下：


从结果看出来几乎完全重合，也就是说明精度提高了很多。