Filter löschen
Filter löschen

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

11 Ansichten (letzte 30 Tage)
guanxu
guanxu am 23 Apr. 2022
Bearbeitet: 埃博拉酱 am 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;%大挠度理论

Antworten (1)

埃博拉酱
埃博拉酱 am 24 Apr. 2022
Bearbeitet: 埃博拉酱 am 24 Apr. 2022
这种简单问题断点调试一下就行了

Kategorien

Mehr zu 常微分方程 finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!