二阶非线性常微分方程组

玉林

我得到的方程组和振动方程所得的方程组类似，
Mq'' +Dq' +Kq=f
其中q是未知向量，而M D K f是q的函数，" ' "表示对时间求导。
因此质量矩阵和力向量都含有非线性项。
能不能用ode 求解？
有没有这方面的示例？

玉林

这个程序该怎么写呢？

qibbxxt

这样的例子很多，把二阶化成两个一阶，就可以了，剩下的帮助文件里面有详细介绍

ChaChing

M, D, K, f是q的函数, 也可以二阶化成两个一阶?
忘干净了, 改天再查下书!:@L

玉林

回复 ChaChing 的帖子

我也一直在考虑这个问题，想遍所有方法了，可是还是没有进展...

玉林

回复 qibbxxt 的帖子

你好，我现在刚刚接触matlab ，需要用matlab 求解非线性微分方程组，但是好几天都没有什么进展，可以帮帮我吗？谢谢，

aspen

参考http://forum.vibunion.com/thread-42367-1-1.html

玉林

回复 aspen 的帖子

我看过那个帖子，但是他的质量矩阵M是常数矩阵，而我的矩阵是未知量的函数，所以我不想对质量矩阵求逆。还是不能解决...

雪缘

玉林 发表于 2010-9-24 22:46

                               
登录/注册后可看大图

回复 aspen 的帖子

我看过那个帖子，但是他的质量矩阵M是常数矩阵，而我的矩阵是未知量的函数，所以我不想 ...

变系数的，自己根据规则推导一下，把二阶微分方程改写为一阶状态方程就可以了

玉林

我编写的程序如下，可以下面的不知该怎么进行了
function [dpdt] = liang2(t,p,dp)
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
syms E I miu   x                               %定义梁常数
syms f                                              %定义载荷
%uu=p(1) vv=p(2)  p(3)=p(3)  q1=p(4)  q2=p(5)   q3=p(6)  q(4)=p(7)
qe=q(4:7);
dqe=dp(4:7);
l=1;
II=[0 -1;1 0];
A=[cos(p(3)) -sin(p(3));sin(p(3)) cos(p(3))];  %转换矩阵
%转换矩阵关于p(3)=theta 的导数
C1=jacobian(A(:,1),p(3));
C2=jacobian(A(:,2),p(3));
C=[C1 C2];
rou=[x 0].';
%定义形函数S
alpha=(2*x-l)/l;
S1=1-3*alpha^2+2*alpha^3;
S2=(alpha-2*alpha^2+alpha^3)*l;
S3=3*alpha^2-2*alpha^3;
S4=(alpha^3-alpha^2)*l;
S=[0 0 0 0;S1 S2 S3 S4];
e=S*qe;
% define stiff matrix
beta=E*I/l;
K=beta*[12 6*l -12 6*l;6*l 4*l^2 -6*l 2*l^2;-12 -6*l 12 -6*l;6*l 2*l^2 -6*l 4*l^2];
kk=zeros(7,7);
kk=sym(kk);
kk(4:7,4:7)=K(:,:);
% define mass matrix
D=eye(2);
B=C*(rou+e);
L=[D B A*S];
P=L.'*L;
M=jifen(P);
%M对时间的导数
dM=zeros(7,7);
dM=sym(dM);
mrrt=zeros(2,2);
mrrt=sym(mrrt);
y1=jifen(S);
mrct=-A*p(3)*y1*qe+C*y1*dqe;
mret=C*dp(3)*y1;
l2=rou.'*S;
y2=jifen(l2);
l3=S.'*S;
mee=jifen(l3);
mcct=2*y2*dqe+2*qe.'*mee*dqe;
l4=S.'*II*S;
y3=jifen(l4);
mcet=dqe.'*y3;
meet=zeros(4,4);
meet=sym(meet);
dM=[mrrt mrct mret;mrct.' mcct mcet;mret.' mcet.' meet];
jaco=sym(zeros(14,7));
dmm=zeros(7,7);
jaco=[-dM dmm];
T=0.5*dp.'*M*dp;                             %动能
qv2=jacobian(T,p);
qv1=dM*dp;
%qv=qv2.'-qv1;                                  %非线性项
%定义载荷
%F=zeros(2),1);
%F=sym(F);
F=[0;(x+l/2)*f/l];
qf1=L.'*F;
qf11=qf1-qv2.';
%方程组
options=odeset('mass',M,'jacobian',jaco);
qf2=M*dp;
qf=sym(zeros(14,1));
qf=[qf11;qf2];
dpdt=inline('')
end

FtpAdmin

上面的是你的模型函数吗？

玉林

回复 FtpAdmin 的帖子

是啊，不过没有写完。已经把质量矩阵和弹性矩阵、力向量都写出来了。但是不知道该怎么写方程...

Happy99

看这麽长的程序着实累, 建议可否贴出方程!?

FtpAdmin

http://forum.vibunion.com/thread-35208-1-1.html

根据帖子中的方法，把方程改写为状态方程的形式，然后就可以求解了

玉林

回复 Happy99 的帖子

方程是根据拉格朗日方程所得。