|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
求解一个微分方程组 Mx''+Cx'+Kx=F
主程序算出了M C K 都为5×5的矩阵, F的表达式随时间t的变量,用ODE求解
主程序
[t,X]=ode45(@(t,X)thermal1(t,X,M,K,C,tt,F1,f1,f2,f3,f4),[0:120],[0 0 0 0 0 0 0 0 0 0]);
F1,f1,f2,f3,f4 为F的表达式
tt = linspace(0,120,120);
dT1=0;
dT2=0;
for ii=1:1:100
dT1=dT1+((-1)^ii/ii^2)*(-ii^2*pi^2*k2/h^2)*exp(-ii^2*pi^2*k2*tt/h^2)*(cos(ii*pi+ii*pi*h1/2/h)-cos(-ii*pi*h1/2/h));
end
MT1=Efc*af2*(w*h1*h/2)*(qs*h/k1)*(-2/pi^2)*dT1;
F1=rho*A*f11*MT1/2/EI;
% test'error'
for ii=1:1:100
dT2=dT2+((-1)^ii/ii^2)*(-ii^2*pi^2*k2/h^2)^2.*exp(-ii^2*pi^2*k2*tt/h^2)*(cos(ii*pi+ii*pi*h1/2/h)-cos(-ii*pi*h1/2/h));
end
MT2=Efc*af2*(w*h1*h/2)*(qs*h/k1)*(-2/pi^2)*dT2;
MTT2=-MT2/2/EI;
MTT1=-MT1/2/EI;
f1=-quadl(@(x)f_fs1(x,1),0,L)*MTT2*rho*A-quadl(@(x)f_fs1(x,1),0,L)*MTT1*Cdamp;
f2=-quadl(@(x)f_fs1(x,2),0,L)*MTT2*rho*A-quadl(@(x)f_fs1(x,2),0,L)*MTT1*Cdamp;
f3=-quadl(@(x)f_fs1(x,3),0,L)*MTT2*rho*A-quadl(@(x)f_fs1(x,3),0,L)*MTT1*Cdamp;
f4=-quadl(@(x)f_fs1(x,4),0,L)*MTT2*rho*A-quadl(@(x)f_fs1(x,4),0,L)*MTT1*Cdamp;
就是在主程序里面先将F求出来了,然后子程序
function dX=thermal1(t,X,M,K,C,tt,F1,f1,f2,f3,f4)
F1= interp1(tt,F1,t);
f1= interp1(tt,f1,t);
f2= interp1(tt,f2,t);
f3= interp1(tt,f3,t);
f4= interp1(tt,f4,t);
dX=[M*X(6:10); F-C*X(6:10)-K*X(1:5)];
这样子求可以不可以 为什么我求出的结果是错误的呢!
F=[F1;f1;f2;f3;f4;];
|
|