在使用ode45函数​求解微分方程出错,无​法执行赋值,因为左侧​和右侧的元素数目不同

조회 수: 2 (최근 30일)
guanxu
guanxu 2022년 4월 23일
편집: 埃博拉酱 2022년 4월 24일
麻烦大家帮忙看看,下面是主程序和子程序
clear
global EI Fx Fy
d=2;
I=pi*d^4/64;
E=200e3;
EI=E*I;
L=1000;
n=100;
Fcr=pi^2*EI/L^2;
Fx=0.99*Fcr;
k1=4*EI/L;
k2=3.5*EI/L;
thita1=pi/5;
% thita20=-thita1*k1/k2;
% Fy=-(k1*thita1+k2*thita20)/L;
% x0=[thita1,k1*thita1/EI];
x1(1)=0;
y1(1)=0;
ds=L/n;
m=10;
for j=1:m
err=0.1;
num=1;
x0=[thita1*j/m, k1*thita1*j/m/EI]
Y(1,:)=x0;
x(1)=0;
Fx=0.99*Fcr;
Fy=0;
while err>0.005 && num<100
x01=x0;
for i=2:n+1
x(i)=x(i-1)+ds;
[Ti,Yi] = ode45(@bend_columnB,[x(i-1) x(i)],x01);
T(i)=Ti(end);
Y(i,:)=Yi(end,:);
x01=Y(i,:);
x1(i)=x1(i-1)+ds*cos(Y(i,1));
y1(i)=y1(i-1)+ds*sin(Y(i,1));
end
err=abs(y1(end));
thita2=Y(end,1);
Fx=Fx+y1(end)*0.025/j*Fcr;
Fy=-(k1*thita1*j/m+k2*thita2)/L;
num=num+1;
end
num_T(j)=num;
dl(j)=L-x1(end);
FFx(j)=Fx;
FFy(j)=Fy;
figure(2)
plot(x1,y1);hold on
end
figure(1)
plotyy(x,Y(:,1),x,Y(:,2));
% figure(2)
% plot(x1,y1);
figure(3)
plot(dl,FFx);hold on
function dy=bend_columnB(x,y)
global k dth_ds20%k=(F/EI)^0.5
dy = zeros(2,1); % a column vector
dy(1)=y(2);
%dy(2)=-k^2*y(1);%小挠度理论
dy(2)=-k^2*sin(y(1))+dth_ds20;%大挠度理论

답변 (1개)

埃博拉酱
埃博拉酱 2022년 4월 24일
편집: 埃博拉酱 2022년 4월 24일
这种简单问题断点调试一下就行了

카테고리

Help CenterFile Exchange에서 App 构建에 대해 자세히 알아보기

제품

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!