实验一 面向微分方程的数值积分法仿真

实验一面向微分方程的数值积分法仿真

一、实验目的

1．掌握数值积分法的基本概念、原理及应用；

2．用龙格-库塔法解算微分方程，增加编写仿真程序的能力； 3．分析数值积分算法的计算步长与计算精度、速度、稳定性的关系； 4. 对数值算法中的“病态问题”进行研究。 二、实验内容

1、已知系统微分方程及初值条件

,(0)1y

t y y =+= 取步长0.1h =，试分别用欧拉方程法和RK4法求2t h =时的y 值，并将求得的值与解析解

()21t y t e t =--比较(将三个解绘于同一坐标中，且用数值进行比较)，说明造成差异的原

因。（①编程完成；②选用MATLAB ode 函数完成。） 程序代码如下：

t0=0; tf=2; h=0.1; y1=1; y2=1; y3=1; t1=0; t2=0; t3=0

n=round(tf-t0)/h; for i=1:n

y1(i+1)=y1(i)+h*(2*h+y1(i)); t1=[t1,t1(i)+h]; end for i=1:n

k1=y2(i)+t2(i);

k2=y2(i)+h*k1/2+t2(i)+h/2; k3=y2(i)+h*k2/2+t2(i)+h/2; k4=y2(i)+h*k3+t2(i)+h;

y2(i+1)=y2(i)+h*(k1+2*k2+2*k3+k4)/6; t2=[t2,t2(i)+h]; end for i=1:n

y3(i+1)=2*exp(t3(i))-t3(i)-1; t3=[t3,t3(i)+h];

end

plot(t1,y1,'r',t2,y2,'g',t3,y3,'k') 实验结果如下；

00.51 1.52 2.5

2

4

6

8

10

12

分析：红线为用欧拉法得到的结果，绿线为用四阶龙格—库塔法得到的结果，蓝线为根据解析方程得到的结果。其差异原因主要有两个：1、二者的方法不同，欧拉法是根据一阶微分方程计算得到的，龙格—库塔法是根据四阶微分方程得到的；2、由于步长取为0.1，所以得到的图像与解析解之间存在差异，若将步长取小，则得到的解将更靠近解析解。 2、已知系统的传递函数为32

40.6

()102722.06

G s s s s =

+++ 在单位阶跃输入下，系统响应的解析解为

1.88 1.88 6.24() 1.84 4.95 1.50.34t t t y t te e e ---=---

试分别用欧拉方程法和RK4法对系统进行仿真（编程完成）：

1）比较两种数值积分解与解析解得逼近程度；（绘图） 程序代码如下：

num=[40.6];den=[1 10 27 22.06]; [a,b,c,d]=tf2ss(num,den);

a=[0 1 0; 0 0 1; -22.06 -27 -10]; b=[0;0;1]; c=[40.6 0 0];

X1=[0;0;0];t=0;Y1=0; X=0; u=1; Y2=0;Y3=0; X2=[0;0;0]; x=0;

h=0.1;

t0=0;

tf=2;

t1=0;t2=0;t3=0;

N=(tf-t0)/h;

for i=1:N

k1=a*X1+b;

k2=b+a*(h*k1/2+X1);

k3=b+a*(h*k2/2+X1);

k4=b+a*(h*k3+X1);

X1=X1+h*(k1+2*k2+2*k3+k4)/6;

Y1=[Y1,c*X1];

t1=[t1,t1(i)+h];

end

for i=1:N

x=X2(:,i)+h*(a*X2(:,i)+b*u);

y=c*x;

X2=[X2,x];

Y2=[Y2,y];

t2=[t2,t2(i)+h];

end

for i=1:N

y=1.84-4.95*i*exp(-1.88*i)-1.5*exp(-1.88*i)-0.34*exp(-6.24*i); Y3=[Y3,y];

t3=[t3,t3(i)+h];

end

plot(t1,Y1,'r',t2,Y2,'g',t3,Y3,'b')

当h=0.01时的结果

0.5

1

1.5

2

2.5

00.20.40.60.811.21.41.61.8

2

当h=0.01时的结果

0.5

1

1.5

2

2.5

分析：这是我得到的结果，发现两个方法得到的结果与实际结果都有较大差距，当是龙格—库塔法更接近实际的结果。

2）改变步长，分析步长对数值解精度的影响；

改变步长后，发现只是两根仿真得到的曲线靠近了，但是与实际曲线仍然是差距很大，这是经过仔细的检查和讨论我觉得程序还是对的，不知道错在哪里了。

3）不断加大步长，分析计算稳定性的变化。 当取h=0.5时，得到的结果：

00.20.40.60.81 1.2 1.4 1.6 1.82

加大步长后结果得到的结果不稳定，不能够很好的对系统进行仿真，另外，由于系统步长选择偏大，根据解析解得到的结果也与实际值有了一定的差距，但是如果步长取得不一样又无法比较。

3、求下图所示系统的阶跃响应()y t 的数值解。（，1v =，1k =，00t =，10f t =0.25h =）分析k 、v 对系统响应的影响。（①编程用RK4求解；②Simulink ）

程序代码如下：

k=1;

a=conv([1 0 0],conv([0.25 1],[0.25 1])) b=[2*k k] X0=[0 0 0 0];

V=1; n=4;T0=0;Tf=10; h=0.25;R=1; b=b/a(1);a=a/a(1);A=a(2: n+1);

A=[rot90(rot90(eye(n-1,n)));-fliplr(A)]; B=[zeros(1,n-1),1]'; m1=length(b);

C=[fliplr(b),zeros(1,n-m1)]; Ab=A-B*C*V; X=X0';y=0;t=T0; N=round(Tf-T0)/h; for i=1:N k1=Ab*X+B*R;

k2=Ab*(X+h*k1/2)+B*R; k3=Ab*(X+h*k2/2)+B*R; k4=Ab*(X+h*k3)+B*R; X=X+h*(k1+2*k2+2*k3+k4)/6;

(r