[已答复] MATLAB ode45求解高阶微分方程组,出现“无法执行赋值,因为左侧和右侧的元素数目不同”问题

[复制链接]
@@@@@BEI 发表于 2021-1-10 11:17:17
本帖最后由 @@@@@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)<-1

    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)<-1

    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)<-1

    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)<-1

    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)<-1

    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])%显示范围设置
























QQ截图20210111103338.png

2 条回复


TouAkira 发表于 2021-1-10 19:56:45
你的自定义函数里面全部参数都缺少赋值,自己没运行过 f_fff( 0, y0 ) 测试能不能正常计算的吗?

@@@@@BEI 发表于 2021-1-11 10:24:37
本帖最后由 @@@@@BEI 于 2021-1-11 10:36 编辑
TouAkira 发表于 2021-1-10 19:56
你的自定义函数里面全部参数都缺少赋值,自己没运行过 f_fff( 0, y0 ) 测试能不能正常计算的吗? ...

微分方程里面本来有时变因素,我把所有参数设为常数后运行程序,程序可以运行,我再研究研究,谢谢
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
热门教程
站长推荐
快速回复 返回顶部 返回列表