欧拉法、RK4、和双线性变换仿真.docx

欧拉法、RK4、和双线性变换仿真实验一、1、欧拉法： 根据欧拉法的计算方法，可以得出待仿真系统的差分方程为： x1(n+1)=x1(n)+h*x2(n); x2(n+1)=x2(n)-h*w*x1(n);所以可得其仿真程序如下：function A err=euler(h) %h=0.01;% w=1; x10=10; x20=0; x1(1)=x10+h*x20; x2(1)=x20-h*w*x10; for i=1:1000 x1(i+1)=x1(i)+h*x2(i); x2(i+1)=x2(i)-h*w*x1(i); end t=0:0.01:10; subplot(2,1,1) plot(t,x1); ylabel(); f=10*cos(t); subplot(2,1,2) plot(t,f); ylabel(); A=x1(101); err=f(101)-x1(101);end输入步长为0.1，0.01和0.001时仿真结果如下：h=0.1h=0.01h=0.0012、RK4法：根据RK4方法编制的仿真程序如下：function y err=rk4_2(h,n) %n=10; %h=0.01; x1_0=10; x2_0=0; k11_0=h*x2_0; k21_0=-h*x1_0; k12_0=h*(x2_0+0.5*k21_0); k22_0=-h*(x1_0+0.5*k11_0); k13_0=h*(x2_0+0.5*k22_0); k23_0=-h*(x1_0+0.5*k12_0); k14_0=h*(x2_0+0.5*k23_0); k24_0=-h*(x1_0+0.5*k13_0); x1(1)=x1_0+(1/6)*(k11_0+2*k12_0+2*k13_0+k14_0); x2(1)=x2_0+(1/6)*(k21_0+2*k22_0+2*k23_0+k24_0); for i=1:100*n k11(i)=h*x2(i); k21(i)=-h*x1(i); k12(i)=h*(x2(i)+0.5*k21(i); k22(i)=-h*(x1(i)+0.5*k11(i); k13(i)=h*(x2(i)+0.5*k22(i); k23(i)=-h*(x1(i)+0.5*k12(i); k14(i)=h*(x2(i)+0.5*k23(i); k24(i)=-h*(x1(i)+0.5*k13(i); x1(i+1)=x1(i)+(1/6)*(k11(i)+2*k12(i)+2*k13(i)+k14(i); x2(i+1)=x2(i)+(1/6)*(k21(i)+2*k22(i)+2*k23(i)+k24(i); end y=x1(100*n+1); t=0:0.01:n; f=10*cos(t); err=f(100*n+1)-y; subplot(2,1,1) plot(t,x1); ylabel(rk4); subplot(2,1,2) plot(t,f); ylabel();end取步长为0.1，0.01，0.001时结果如下：h=0.1h=0.01h=0.001通过对以上两种算法的结果的观察与比较我们发现：首先，对于零阶稳定的原系统，由于计算机的舍入误差和算法的阶段误差的存在，不可能得到一个稳定的数值解；其次，当h较小时，RK4的方法的效果更好一些。实验二、1、 双线性替换法：双线性替换法的仿真程序如下：function A=dlinear(t,k) %t=0.1; %k=1; x0=0; x01=0; for j=1:100 x(j)=1; end y0=0; y01=0; y(1)=(t2)*x(1)+2*(t2)*x0+(t2)*x01-2*(t*k)2-4)*y0-(t*k-2)2*y01)/(t*k+2)2); y(2)=(t2)*x(2)+2*(t2)*x(1)+(t2)*x0-2*(t*k)2-4)*y(1)-(t*k-2)2*y0)/(t*k+2)2); for i=3:100 y(i)=(t2)*x(i)+2*(t2)*x(i-1)+(t2)*x(i-2)-2*(t*k)2-4)*y(i-1)-(t*k-2)2*y(i-2)/(t*k+2)2); end h=0.1:0.1:10; subplot(2,1,1) plot(h,y); ylabel(); s=tf(s); G=1/(s+k)2; subplot(2,1,2); step(G); ylabel();end当K=1时，取步长分别为1，0.1，0.01，结果如下：步长为1步长为0.1步长为0.01当K=10时，取步长为1，0.01，0.001结果如下：步长为1步长为0.01步长为0.0012、 根匹配法：其程序设计如下：function A=rootmatch(t,k) %t=0.1; %k=1; xend=1/k2; m=exp(-k*t); x0=0; x01=0; for j=1:100 x(j)=1; end y0=0; y01=0; y(1)=xend*(1-m)2)*x01+2*m*y0-(m2)*y01; y(2)=xend*(1-m)2)*x0+2*m*y(1)-(m2)*y0; for i=3:100 y(i)=xend*(1-m)2)*x(i-2)+2*m*y(i-1)-(m2)*y(i-2); end h=0.01:0.01:1; subplot(2,1,1) plot(h,y); ylabel(); s=tf(s); G=1/(s+k)2; subplot(2,1,2); step(G); ylabel(S)end当k=1时，取步长分别为