欧拉法、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.1

h=0.01

h=0.001

2、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('½âÎö½â');

end

取步长为0.1，0.01，0.001时结果如下：

h=0.1

h=0.01

h=0.001

通过对以上两种算法的结果的观察与比较我们发现：首先，对于零阶稳定的原系统，由于计算机的舍入误差和算法的阶段误差的存在，不可能得到一个稳定的数值解；其次，当h较小时，RK—4的方法的效果更好一些。

实验二、

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)=((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('ÊýÖµ½â');

 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.001

2、根匹配法：

其程序设计如下：

function A=rootmatch(t,k)

  %t=0.1;

  %k=1;

  xend=1/k^2;

  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-(m^2)*y01;

  y(2)=xend*((1-m)^2)*x0+2*m*y(1)-(m^2)*y0;

  for i=3:100

    y(i)=xend*((1-m)^2)*x(i-2)+2*m*y(i-1)-(m^2)*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时，取步长分别为1，0.1，0.01，得结果如下图：

步长为1

步长为0.1

步长为0.01

当K=10时，取步长分别为0.1，0.01，0.001，得结果如下：

步长为0.1

步长为0.01

步长为0.001

  通过对结果的观察，我们发现K对系统的稳定性和仿真步长有着很大的影响。