oneonly 发表于 2015-10-13 16:26
这就是那本结构动力学，我的方程以化成书上的这种形式，
k1=h*(C*y(:,n)+F);
k2=h*(C*(y(:,n)+ ...
这是标准的定不长4阶龙哥库塔法，从效果上来说，其未必会比ode45好
从算法本身来说ode45就有了极大的改善，而且相关的代码是matlab公司多次优化过的
至于调整问题我觉得看看两个方面吧，一个是程序本身有没有问题
我很早以前发过相关的代码，你看看http://forum.vibunion.com/thread-17615-1-1.html
另一方面看看输出的激励变化规律看看
oneonly 发表于 2015-10-13 20:29
happy老师，晚上我出了一组数据y（107，：）=0 -1129183048.93119 1.48637526202825e+32 ...
这个应该是结果发散了
oneonly 发表于 2015-10-9 10:54
额。。。以前师兄是用newmark方法做的，现在呢 是换一种算法，然后在进行后续的分析。现在单单是解这一个 ...
http://forum.vibunion.com/thread-138097-1-1.html
happy 发表于 2015-10-14 09:50
这个应该是结果发散了
以前我也考虑用matlab内部的ode45，但是我的是个变系数矩阵的（K C M F）都是随时间变化的，项数也比较多（感觉应该用循环调用）不怎么会编写程序，所以就自己编的龙格库塔的程序，有一点疑惑就是 （例如）k2=h*f(yn+1/2*k1，tn+1/2*h)，这里面有个时间差tn不知道怎么处理(时间变化，我的系数也变了)，
最后，我就假定时间t=0.3求出 K C M F，并且然他们为定值，运用龙格库塔方法计算，出来的数就非常大（都200多次方），并且后来全是NAN。。。
oneonly 发表于 2015-10-14 17:28
以前我也考虑用matlab内部的ode45，但是我的是个变系数矩阵的（K C M F）都是随时间变化的，项数也比较多 ...
tn是离散点的时间点，也就是yn所对应的时间点
出来NAN意味这出现了0/0,0*无穷，等运算
mxlzhenzhu 发表于 2015-10-14 11:46
http://forum.vibunion.com/thread-138097-1-1.html
这是我的变系数动力学方程：http://forum.vibunion.com/thread-138097-1-1.html按照书中308页下边的那种龙格库塔的方法，结合您的帖子上面的程序，不做变换化为一阶方程，而直接用龙哥-库塔对其求解。我的编程思路大致如下： clcclear %%%%%%%%%%%%%%%%%%%%%%% Ixy=0.81*10^-8; Ixz=0.56*10^-2;Iyx=0.81*10^-8; Iyy=0.66432; Iyz=1.166*10^-7;Izx=0.56*10^-2; Izy=1.166*10^-7; Izz=0.64; 密度，质量 转动惯量等已知是的赋值 %%%%%%%%%%%%%%%%%%%%%%%%%%%% tt=0:0.01:1; %%%%%%这里也试过0.001,0.0001等for ij=1:length(tt) t=tt(ij)dt=0.01;%%%%%%%%%步长也试过0.1,0.0001等 。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。 M= K= F= C= 。。。。。。。。。。。。。。。。。。。。。。 。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。%%%%%句号中间是求M K F C （都是关于t的函数）Inv_M=inv(M); if ij==1q=zeros(112,1);%%%%初始位移q1d=zeros(112,1); %%%%初始速度q2d=Inv_M*(QA-K*q-C*q1d); %%%%初始加速度 endif ij>1qq=q;qq1d=q1d;qq2d=q2d; Kd1=Inv_M*(-K*qq-C*qq1d+F);%Kd2=Inv_M*(-K*(qq+0.5*dt*qq1d)-C*(qq1d+0.5*dt*Kd1)+F);Kd3=Inv_M*(-K*(qq+0.5*dt*qq1d+0.25*dt^2*Kd1)-C*(qq1d+0.5*dt*Kd2)+F);Kd4=Inv_M*(-K*(qq+dt*qq1d+0.5*dt^2*Kd2)-C*(qq1d+dt*Kd3)+F); update_disp=qq+dt*qq1d+(dt^2)/6*(Kd1+Kd2+Kd3);%update_velo=qq1d+dt/6*(Kd1+2*Kd2+2*Kd3+Kd4);update_acc =Kd1; q=update_disp; q1d=update_velo;q2d=update_acc; q107=q(107); %%%x位移 q108=q(108); %%%y位移 q109=q(109); %%%y位移 qqq1(ij)=q107; qqq2(ij)=q108; qqq3(ij)=q109; end end figure(1)plot(tt,qqq1,'b')%%%x位移 title('λòÆ')xlabel('ê±¼ät/s')ylabel('λòÆ/m') figure(2)plot(tt,qqq2,'m')%%%y位移title('Ëù¶è')xlabel('ê±¼ät/s')ylabel('λòÆ/m') figure(3)plot(tt,qqq3,'m')%%%z位移title('Ëù¶è')xlabel('ê±¼ät/s')ylabel('λòÆ/m')
本帖最后由 oneonly 于 2015-10-18 20:08 编辑
happy 发表于 2015-10-16 08:21
tn是离散点的时间点，也就是yn所对应的时间点
出来NAN意味这出现了0/0,0*无穷，等运算
happy 老师您好有时间能帮看一下这个文档吗，谢谢谢谢
本帖最后由 oneonly 于 2015-10-18 20:07 编辑
mxlzhenzhu 发表于 2015-10-14 11:46
http://forum.vibunion.com/thread-138097-1-1.html
老师您好，那边老实需要审核，我就在这里给你回复吧，希望能看到，谢谢谢谢谢谢
oneonly 发表于 2015-10-18 19:48
这是我的变系数动力学方程：http://forum.vibunion.com/thread-138097-1-1.html按照书中308页下边的那 ...
你用该方法求解已经很长时间了，相信代码本身方面为题应该不大
如果你真要用这个方法实现求解，给你个建议吧
去找一本微粉方程数值方法方面的书，有深入将RK法的那种
好好研究一下里边关于收敛性和稳定性的内容
然后在看看你的模型是否在RK法的绝对稳定域内。如果不在这个范围内的，估计是无法收敛或稳定的
那么就只能更换算法或者针对你的具体问题改进RK法