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

28 views (last 30 days)
guanxu on 23 Apr 2022
Edited: 埃博拉酱 on 24 Apr 2022

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;%大挠度理论

Edited: 埃博拉酱 on 24 Apr 2022

### Categories

Find more on 数据定义基础知识 in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!