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

欧拉法、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;%2?3¤ 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???a); A=x1(101); err=f(101)-x1(101);end输入步长为0.1，0.01和0.001时仿真结果如下：h=0.1h=0.01h=0.0012、RK—4法：根据RK—4方法编制的仿真程序如下：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(?a???a);end取步长为0.1，0.01，0.001时结果如下：h=0.1h=0.01h=0.001通过对以上两种算法的结果的观察与比较我们发现：首先，对于零阶稳定的原系统，由于计算机的舍入误差和算法的阶段误差的存在，不可能得到一个稳定的数值解；其次，当h较小时，RK—4的方法的效果更好一些。实验二、双线性替换法：双线性替换法的仿真程序如下：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)=((t^2)*x(1)+2*(t^2)*x0+(t^2)*x01-2*((t*k)^2-4)*y0-(t*k-2)^2*y01)/((t*k+2)^2); y(2)=((t^2)*x(2)+2*(t^2)*x(1)+(t^2)*x0-2*((t*k)^2-4)*y(1)-(t*k-2)^2*y0)/((t*k+2)^2); for i=3:100 y(i)=((t^2)*x(i)+2*(t^2)*x(i-1)+(t^2)*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(êy?μ?a); s=tf(s); G=1/(s+k)^2; subplot(2,1,2); step(G); ylabel(?a???a);end当K=1时，取步长分别为1，0.1，0.01，结果如