yuedu baben.txt

基于matlab开发软件 的选择输入类型的数值计算仿真程序 迭代方法包括 1可自启动的4阶龙格库塔法 2需要使用其他方法启动的阿达姆斯法

function zjmain
clear                                   %清除工作空间
clc                                     %清除命令窗口
clf
zj1=input('请选择:0(输入传递函数求解状态空间方程)1(输入状态空间方程系数)');
if(zj1==0)
    fm=input('请输入分母系数:');            %人机交互输入分母系数矩阵
    fz=input('请输入分子系数:');            %人机交互输入分子系数矩阵
    n=length(fm)-1;                         %求n
    m=length(fz)-1;                         %求m
    a=zeros(n,n);                           %初始化a
    for i=1:n                               %求a
        if(i>1)
            a(i-1,i)=1;
        end
        a(n,i)=-fm(n-i+2);
    end
    b=zeros(n,1);b(n,1)=1;                  %求b
    c=zeros(1,n);                           %初始化c     
    if (m<n)                                % m<n 情况下求c,d 
        for i=1:m+1
            c(1,i)=fz(m-i+2);
        end
        d=0;
    elseif (m==n)                          % m=n 情况下求c,d
        for i=1:m
            c(1,i)=fz(m-i+2)-fz(1).*fm(n-i+2);
        end
        d=fz(1);
    end
    hui2=input('系统传递函数为：0(开环传函),1(闭环传函):');
    if(hui2==0)
        v=input('请输入反馈系数');
        a1=a'-(c')*v*(b');b1=c';c1=b';d1=d;
    elseif(hui2==1)
        a1=a';b1=c';c1=b';d1=d;
    end   
    fprintf('能观型状态空间方程组：\n');
    for i=1:n
        if (i==n)
            fprintf('x''  = |');
        else
            fprintf('      |');
        end
        fprintf('%5d',a1(i,:));
        if (i==n)
            fprintf('|x + |%5d|u\n\n',b1(i,1));
        else
            fprintf('|    |%5d| \n\n',b1(i,1));
        end
    end
    fprintf('y = [');
    fprintf('%5d',c1);
    if (d1~=0)
        fprintf(']x + %3du\n\n',d1);
    else
        fprintf(']x\n\n');
    end
elseif (zj1==1)
    a11=input('a=');
    b11=input('b=');
    c11=input('c=');
    d11=input('d=');
    hui2=input('系统传递函数为：0(开环传函),1(闭环传函):');
    if(hui2==0)
        v=input('请输入反馈系数');
        a1=a11-b11*v*c11;b1=b11;c1=c11;d1=d11;
    elseif(hui2==1)
        a1=a11;b1=b11;c1=c11;d1=d11;
    end
end
zj2=input('已知初值是:0(y,u),1(x,u),2(x)?');
if (zj2==0)
    y0=input('请输入y的初值:');                                  %输入y的初值     
    u00=input('请输入u的初值:');                                  %输入u的初值
    x0=zeros(1,n);
    for i=n:-1:1                                                 %计算初值
        x0(i)=y0(i);
        for j=i:n-1
            if(m==0)
                x0(i)=x0(i)-a1(n+j-2,n).*y0(n-j+i);
            else
                x0(i)=x0(i)-a1(n+j-2,n).*y0(n-j+i)-b1(m+j-1).*u00(m-j+i);
            end
        end
    end
    u0=u00(m);
    fprintf('初值为:');x0=x0';x0
elseif (zj2==1)
    u0=input('请输入u的初值:');
    x0=input('请输入x的初值:');
elseif (zj2==2)
    x0=input('请输入x的初值:');
    u0=0;
end
zj4=input('请选择数值求解方法:0(单步长运算),1(不同步长运算)');
if(zj4==0)
    fprintf('四阶龙格库塔法,阿达姆斯法两种方法比较:\n');
    h=input('请输入仿真步长: ');
    n=input('请输入迭代次数: ');
    [x1,y1]=zjode4(a1,b1,c1,d1,u0,h,n,x0);
    [x2,y2]=zjadms(a1,b1,c1,d1,u0,h,n,x0);
    fprintf('四阶龙格库塔法初值加前10步的迭代结果(第一段状态变量，第二段输出)\n');
    x1(:,1:11)
    fprintf('y=\n    ');
    fprintf('%.4f    ',y1(:,1:11));
    fprintf('\n阿达姆斯法初值加前10步的迭代结果(第一段状态变量，第二段输出)\n');
    x2(:,1:11)
    fprintf('y=\n    ');
    fprintf('%.4f    ',y2(:,1:11));
    t=0:h:h*n;
    subplot(2,2,1);
    plot(t,y1,'m')
    title('四阶龙格库塔法输出结果)');
    xlabel('时间t');ylabel('输出y');
    subplot(2,2,2)
    plot(t,y2,'r')
    title('阿达姆斯法输出结果');
    xlabel('时间t');ylabel('输出y');
    subplot(2,2,3);
    plot(t,x1,'b')
    title('四阶龙格库塔法状态变量');
    xlabel('时间t');ylabel('状态变量x');
    subplot(2,2,4)
    plot(t,x2,'k')
    title('阿达姆斯法状态变量');
    xlabel('时间t');ylabel('状态变量x');
elseif(zj4==1)
    h1=input('请输入仿真步长数组: '); 
    n=input('请输入迭代次数: ');
    for i=1:length(h1)
        h=h1(i);
        [xx1,yy1]=zjode4(a1,b1,c1,d1,u0,h,n,x0);
        [xx2,yy2]=zjadms(a1,b1,c1,d1,u0,h,n,x0);
        t=0:h:h*n;
        figure(01)
        if(rem(i,2)==1)
            plot(t,yy1,'b')
        else
            plot(t,yy1,'--r')
        end
        xlabel('时间t');ylabel('输出y');
        title('四阶龙格库塔法输出结果');
        hold on
        figure(02)
        if(rem(i,2)==1)
            plot(t,yy2,'b')
        else
            plot(t,yy2,'r')
        end
        xlabel('时间t');ylabel('输出y');
        title('阿达姆斯法输出结果');
        hold on
        figure(03)
        if(rem(i,2)==1)
            plot(t,xx1,'ob')
        else
            plot(t,xx1,'--r')
        end
        xlabel('时间t');ylabel('状态变量x');
        title('四阶龙格库塔法状态变量');
        hold on
        figure(04)
        if(rem(i,2)==1)
            plot(t,xx2,'b')
        else
            plot(t,xx2,'--r')
        end
        xlabel('时间t');ylabel('状态变量x');
        title('阿达姆斯法状态变量');
        hold on
        fprintf('\n仿真步长 %f\n四阶龙格库塔法初值加前10步的迭代结果(第一段状态变量，第二段输出)\n',h);
        xx1(:,1:11)
        fprintf('y=\n    ');
        fprintf('%.4f    ',yy1(:,1:11));
        fprintf('\n阿达姆斯法初值加前10步的迭代结果(第一段状态变量，第二段输出)\n');
        xx2(:,1:11)
        fprintf('y=\n    ');
        fprintf('%.4f    ',yy2(:,1:11));
    end
end
k=input('\n是否重新输入:1(是) 0(否)?');
if(k==1)
    zjmain                                                          %重新执行函数
end






function [x1,y1]=zjode4(a,b,c,d,u,h,n,x0)
x=x0;
y=c*x;
for i=1:n
    k1=a*x(:,i)+b*u;
    k2=(a*(x(:,i)+0.5*h*k1)+b*u);
    k3=(a*(x(:,i)+0.5*h*k2)+b*u);
    k4=(a*(x(:,i)+h*k3)+b*u);
    x(:,i+1)=x(:,i)+h/6*(k1+2*k2+2*k3+k4);
    y(i+1)=c*x(:,i+1);
end
x1=x;
y1=y;





function [x2,y2]=zjadms(a,b,c,d,u,h,n,x0)
x=x0;
y=c*x;
for i=1:n
    if (i==1)
        x(:,i+1)=x(:,i)+h*(a*x(:,i)+b*u);
    else
        x(:,i+1)=x(:,i)+h/2*(3*(a*x(:,i)+b*u)+(a*x(:,i-1)+b*u));
    end
    y(i+1)=c*x(:,i+1);
end
x2=x;
y2=y;