Question

MATLAB 程序求助

%define the sequence of the program
[A,B,U,Y,M]=define_matrix();
[H,D,I,R1,C2,C1,R2,P,Cd,RO,Q0,X1,P0,G]=input_parameter();
[X2,C,A1]=calculation_parameter(Q0,Cd,P,R2,D,RO,G);
[M]=calculation_matrix(M,A,B,U,Y,I,R1,C2,C1,R2,P0,P,X2,X2,A1,Q0,H,C);
graph_plot(M);

function [A,B,U,Y,M]=define_matrix()
%A_state_matrix
A=zeros(3,3);
%B=input_matrix
B=zeros(3,2);
%U_input_vector
U=zeros(2,1);
%Y_state_matrix
Y=zeros(3,1);
%M_output_matrix
M=zeros(1000,4);

function [H,D,I,R1,C2,C1,R2,P,Cd,RO,Q0,X1,P0,G]=input_parameter()
%步长H(s)，阀芯直径D（cm）,惯性质量I(kg)，阀孔阻尼R1（bar_s/cm3）
%弹簧柔度C2(cm/kg），管道液容C1（cm3/bar），泄油液阻R2(bar_s/cm3)
%设定压力P（bar），初始压力P0（bar），流量系数Cd，油液密度RO（kg/cm3）
%泵理论流量Q0（cm3/s），阀口遮盖量X1（cm），重力加速度G（cm/s2）
H=6e-5; D=1.2; I=6.3e-5; R1=0.4; R2=1.5; C2=0.02;
C1=0.08; P=30; Cd=0.7; RO=0.0009; Q0=460; X1=0.14; P0=6; G=981;

function [X2,C,A1]=calculation_parameter(Q0,Cd,P,R2,D,RO,G)
PI=3.1415926;
X2=(Q0-P/R2)/(Cd*PI*D*sqrt(2*G*P/RO));
C=-Cd*PI*D*sqrt(2*G/RO);
A1=PI*D*D/4;

function [M]=calculation_matrix(M,A,B,U,Y,I,R1,C2,C1,R2,P0,P,X1,
X2,A1,Q0;H,C)
T=0;
A(1,1)=-A1*A1*R1/I; A(1,2)=-1/C2; A(1,3)=A1/C1;
A(2,1)=1/I; A(3,1)=-A1/I; A(3,3)=-1/(R2*C1);
B(1,1)=-1; B(3,2)=1;
U(1)=P*A1-(X1+X2)/C2; U(2)=Q0;
Y(1)=0; Y(2)=0; Y(3)=P0*C1;
M(1,1)=0; M(1,2)=Y(1)/I; M(1,3)=Y(2); M(1,4)=P0;
for j=1:1000,
[Y,T]=R_T(A,B,U,Y,H,X1,C,C1,T);
M(j+1,1)=T; M(j+1,2)=Y(1)/I; M(j+1,3)=Y(2); M(j+1,4)=Y(3)/C1;
end

function [Y,T]=R_T(A,B,U,Y,H,X1,C,C1,T)
%constrain
if Y(2)<0
Y(2)=0;
end
if Y(2)==0&&Y(1)<0
Y(1)=0;
end

%nonlinear_situation
if Y(2)>X1
A(3,3)=A(3,3)+C*(Y(2)-X1)*sqrt(Y(3)/C1);
end

W=1/3*[0.5,1,1,0.5];
Q=[0.5,0.5,1,0];
TW=T; BB=Y; YW=Y;
for i=1:4
Dy=A*Y+B*U;
T=TW+W(i)*H;
Y=YW+H*Q(i).*Dy;
BB=BB+H*W(i).Dy;
end
Y=BB;
if Y(2)<0
Y(2)=0;
else
if Y(2)==0&&Y(1)<0
Y(1)=0;
end
end

function graph_plot(M)
%%plot_graph
figure(1);
plot(M(:,1),M(:,2));
figure(2);
plot(M(:,1),M(:,3));
figure(3);
plot(M(:,1),M(:,4));

Answer 1

??? Error: File: untitled.m Line: 8 Column: 1
Function definitions are not permitted at the prompt or in scripts.

这个是提示说你写的函数文件要单独建立在一个函数文件中，而不能直接将函数放在脚本文件中。具体方法如下：
一、先建立主程序 main.m
%define the sequence of the program
[A,B,U,Y,M]=define_matrix();
[H,D,I,R1,C2,C1,R2,P,Cd,RO,Q0,X1,P0,G]=input_parameter();
[X2,C,A1]=calculation_parameter(Q0,Cd,P,R2,D,RO,G);
[M]=calculation_matrix(M,A,B,U,Y,I,R1,C2,C1,R2,P0,P,X2,X2,A1,Q0,H,C);
graph_plot(M);

二、建立函数文件1 define_matrix.m
function [A,B,U,Y,M]=define_matrix()
%A_state_matrix
A=zeros(3,3);
%B=input_matrix
B=zeros(3,2);
%U_input_vector
U=zeros(2,1);
%Y_state_matrix
Y=zeros(3,1);
%M_output_matrix
M=zeros(1000,4);

三、建立函数文件2 input_parameter.m
function [H,D,I,R1,C2,C1,R2,P,Cd,RO,Q0,X1,P0,G]=input_parameter()
%步长H(s)，阀芯直径D（cm）,惯性质量I(kg)，阀孔阻尼R1（bar_s/cm3）
%弹簧柔度C2(cm/kg），管道液容C1（cm3/bar），泄油液阻R2(bar_s/cm3)
%设定压力P（bar），初始压力P0（bar），流量系数Cd，油液密度RO（kg/cm3）
%泵理论流量Q0（cm3/s），阀口遮盖量X1（cm），重力加速度G（cm/s2）
H=6e-5; D=1.2; I=6.3e-5; R1=0.4; R2=1.5; C2=0.02;
C1=0.08; P=30; Cd=0.7; RO=0.0009; Q0=460; X1=0.14; P0=6; G=981;

四、建立函数文件3 calculation_parameter.m
function [X2,C,A1]=calculation_parameter(Q0,Cd,P,R2,D,RO,G)
PI=3.1415926;
X2=(Q0-P/R2)/(Cd*PI*D*sqrt(2*G*P/RO));
C=-Cd*PI*D*sqrt(2*G/RO);
A1=PI*D*D/4;

五、建立函数文件4 calculation_matrix.m
function [M]=calculation_matrix(M,A,B,U,Y,I,R1,C2,C1,R2,P0,P,X1,
X2,A1,Q0;H,C)
T=0;
A(1,1)=-A1*A1*R1/I; A(1,2)=-1/C2; A(1,3)=A1/C1;
A(2,1)=1/I; A(3,1)=-A1/I; A(3,3)=-1/(R2*C1);
B(1,1)=-1; B(3,2)=1;
U(1)=P*A1-(X1+X2)/C2; U(2)=Q0;
Y(1)=0; Y(2)=0; Y(3)=P0*C1;
M(1,1)=0; M(1,2)=Y(1)/I; M(1,3)=Y(2); M(1,4)=P0;
for j=1:1000,
[Y,T]=R_T(A,B,U,Y,H,X1,C,C1,T);
M(j+1,1)=T; M(j+1,2)=Y(1)/I; M(j+1,3)=Y(2); M(j+1,4)=Y(3)/C1;
end

六、建立函数文件5 R_T.m
function [Y,T]=R_T(A,B,U,Y,H,X1,C,C1,T)
%constrain
if Y(2)<0
Y(2)=0;
end
if Y(2)==0&&Y(1)<0
Y(1)=0;
end

%nonlinear_situation
if Y(2)>X1
A(3,3)=A(3,3)+C*(Y(2)-X1)*sqrt(Y(3)/C1);
end

W=1/3*[0.5,1,1,0.5];
Q=[0.5,0.5,1,0];
TW=T; BB=Y; YW=Y;
for i=1:4
Dy=A*Y+B*U;
T=TW+W(i)*H;
Y=YW+H*Q(i).*Dy;
BB=BB+H*W(i).Dy;
end
Y=BB;
if Y(2)<0
Y(2)=0;
else
if Y(2)==0&&Y(1)<0
Y(1)=0;
end
end

七、建立函数文件6 graph_plot.m
function graph_plot(M)
%%plot_graph
figure(1);
plot(M(:,1),M(:,2));
figure(2);
plot(M(:,1),M(:,3));
figure(3);
plot(M(:,1),M(:,4));

追问

太感谢了，是否可以加我为好友，以后有问题了可以方便请教你

追答

可以买一本MATLAB的教材自学，MATLAB不难~天天向上~

Answer 2

Answer 3

klop>FPHI
at
26
CK=CK*(1.0-Z^2);
改成：
CK=CK.*(1.0-Z.^2);
.运算是matlab特有的，可以查一下用法，也是初学者经常犯的错误
z^2是直接平方是矩阵平方，而只有方矩阵才能平方，而z.^2是对数组或矩阵每个元素分别平方，跟矩阵平方不一样
同样a*b是矩阵乘，a.*b是数组或矩阵的对应元素相乘

Answer 4

yn(i,a)<0.5
就是yn矩阵第i行第a列元素小于0.5

追问

??? Error: File: untitled.m Line: 8 Column: 1
Function definitions are not permitted at the prompt or in scripts.
这个是怎么回事