本帖最后由 @@@@@BEI 于 2021-1-11 10:38 编辑

下面是程序的主要部分，方程本来有间隙函数、时变刚度、变外力等等，为了运行程序，参数暂时全部取常数，再试运行程序。就是微分方程总是出现这个问题，还有问题就是：“输入向量与输出向量不一样”，或者“返回的向量和初始条件向量的元素数目必须相同”，方程本来是15阶的二阶微分方程，降阶后有30个一元方程组，用ode45龙格库塔计算，就出现维度不一致的问题·，搞了大半天也不明白哪里出了问题，求教求教。

%函数定义

% function dy=f_ffff(t,y)

dy=zeros(30,1);

z_3=f_3(y);

z_5=f_5(y);

z_7=f_7(y);

z_21=f_21(y);

z_23=f_23(y);

global k11 k12 k13 k14 k21 k22 k25 k28 k213 k31 k33 k36 k39 k314 k41 k44 k47 k410 k415

global k52 k55 k511 k63 k66 k611 k74 k77 k711 k82 k88 k812 k93 k99 k912 k104 k1010 k1012

global k115 k116 k117 k1111 k128 k129 k1210 k1212 k132 k1313 k143 k1414 k154 k1515

global c11 c12 c13 c14 c21 c22 c25 c28 c213 c31 c33 c36 c39 c314 c41 c44 c47 c410 c415

global c52 c55 c511 c63 c66 c611 c74 c77 c711 c82 c88 c812 c93 c99 c912 c104 c1010 c1012

global c115 c116 c117 c1111 c128 c129 c1210 c1212 c132 c1313 c143 c1414 c154 c1515

global F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 F13 F14 F15

%无量纲刚度矩阵

k11=1;

k12=1;

k13=1;

k14=1;

k21=1;

k22=1;

k25=1;

k28=1;

k213=1;

k31=1;

k33=1;

k36=1;

k39=1;

k314=1;

k41=1;

k44=1;

k47=1;

k410=1;

k415=1;

k52=1;

k55=1;

k511=1;

k63=1;

k66=1;

k611=1;

k74=1;

k77=1;

k711=1;

k82=1;

k88=1;

k812=1;

k93=1;

k99=1;

k912=1;

k104=1;

k1010=1;

k1012=1;

k115=1;

k116=1;

k117=1;

k1111=1;

k128=1;

k129=1;

k1210=1;

k1212=1;

k132=1;

k1313=1;

k143=1;

k1414=1;

k154=1;

k1515=1;

%无量纲阻尼矩阵

c11=1;

c12=1;

c13=1;

c14=1;

c21=1;

c22=1;

c25=1;

c28=1;

c213=1;

c31=1;

c33=1;

c36=1;

c39=1;

c314=1;

c41=1;

c44=1;

c47=1;

c410=1;

c415=1;

c52=1;

c55=1;

c511=1;

c63=1;

c66=1;

c611=1;

c74=1;

c77=1;

c711=1;

c82=1;

c88=1;

c812=1;

c93=1;

c99=1;

c912=1;

c104=1;

c1010=1;

c1012=1;

c115=1;

c116=1;

c117=1;

c1111=1;

c128=1;

c129=1;

c1210=1;

c1212=1;

c132=1;

c1313=1;

c143=1;

c1414=1;

c154=1;

c1515=1;

%外力

F1=200;

F2=200;

F3=200;

F4=200;

F5=0;

F6=0;

F7=0;

F8=0;

F9=0;

F10=0;

F11=100;

F12=100;

F13=-100;

F14=-100;

F15=-100;

%系统偏微分方程

dy(1)=y(2);

dy(2)=-k11*y(1)-k12*z_3-k13*z_5-k14*z_7-c11*y(2)-c12*y(4)-c13*y(6)-c14*y(8)+F1;

dy(3)=y(4);

dy(4)=-k21*y(1)-c21*y(2)-k22*z_3-c22*y(4)-k25*y(9)-c25*y(10)-k28*y(9)-c28*y(10)-k213*y(25)-c213*y(26)+F2;

dy(5)=y(6);

dy(6)=-k31*y(1)-c31*y(2)-k33*z_5-c33*y(6)-k36*y(11)-c36*y(12)-k39*y(17)-c39*y(18)-k314*y(27)-c314*y(28)+F3;

dy(7)=y(8);

dy(8)=-k41*y(1)-c41*y(2)-k44*z_7-c44*y(8)-k47*y(13)-c47*y(14)-k410*y(19)-c410*y(20)-k415*y(27)-c415*y(28)+F4;

dy(9)=y(10);

dy(10)=-k52*z_3-c52*y(4)-k55*y(9)-c55*y(10)-k511*z_21-c511*y(22)+F5;

dy(11)=y(12);

dy(12)=-k63*z_5-c63*y(6)-k66*y(11)-c66*y(11)-k611*z_21-c611*y(22)+F6;

dy(13)=y(14);

dy(14)=-k74*z_7-c74*y(8)-k77*y(13)-c77*y(14)-k711*z_21-c711*y(22)+F7;

dy(15)=y(16);

dy(16)=-k82*z_3-c82*y(4)-k88*y(15)-c88*y(16)-k812*z_23-c812*y(24)+F8;

dy(17)=y(18);

dy(18)=-k93*z_5-c52*y(6)-k55*y(17)-c55*y(18)-k511*z_23-c511*y(24)+F9;

dy(19)=y(20);

dy(20)=-k104*z_7-c104*y(8)-k1010*y(19)-c1010*y(20)-k1012*z_23-c1012*y(24)+F10;

dy(21)=y(22);

dy(22)=-k115*y(9)-c115*y(10)-k116*y(11)-c116*y(12)-k117*y(13)-c117*y(14)-k1111*z_21-c1111*y(22)+F11;

dy(23)=y(24);

dy(24)=-k128*y(15)-c128*y(16)-k129*y(17)-c129*y(18)-k1210*y(19)-c1210*y(20)-k1212*z_23-c1212*y(24)+F12;

dy(25)=y(26);

dy(26)=-k132*z_3-c132*y(4)-k1313*y(25)-c1313*y(26)+F13;

dy(27)=y(28);

dy(28)=-k143*z_5-c143*y(6)-k1414*y(27)-c1414*y(28)+F14;

dy(29)=y(30);

dy(30)=-k154*z_7-c132*y(8)-k1515*y(29)-c1515*y(30)+F15;

%间隙函数

function fq_3 =f_3(y)

if y(3)>1

fq_3=y(3)-1;

elseif y(3)

fq_3=y(3)+1;

else

fq_3=0;

end

function fq_5 =f_5(y)

if y(5)>1

fq_5=y(5)-1;

elseif y(5)

fq_5=y(5)+1;

else

fq_5=0;

end

function fq_7 =f_7(y)

if y(7)>1

fq_7=y(7)-1;

elseif y(7)

fq_7=y(7)+1;

else

fq_7=0;

end

function fq_21 =f_21(y)

if y(21)>1

fq_21=y(21)-1;

elseif y(21)

fq_21=y(21)+1;

else

fq_21=0;

end

function fq_23 =f_23(y)

if y(23)>1

fq_23=y(23)-1;

elseif y(23)

fq_23=y(23)+1;

else

fq_23=0;

end

用户程序：

T=2*pi;

tspan=0:T/100:200*T;

y0=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

options=odeset('RelTol',10^-3,'AbsTol',10^-3);

[t,y]=ode45(@f_fffff,tspan,y0);

plot(t,y(:,1));%绘制 y 的时间历程

xlabel('t')%横轴为 t

ylabel('y')%纵轴为 y

grid;%显示网格线

axis([460 500 -Inf Inf])%显示范围设置

2021-1-11 10:37 上传