高阶系数矩阵的二阶微分方程

inoran

各位高手好 ：
小弟是MATLAB初学
看了HELP 没有说到如何用矩阵的写法解
二阶微分方程式,所以我自己修改了书本上面一阶的写法

例如：  [M]y''+[K]y=[f11]*sin(w*t)   [M]与[K]如果是6*6  [f11]6*1

这要我要写成程式的话我知道要先降皆,以下我先假设sin(w*t)是1  
副程式如下：

function xdot=inoran(t,x)
global M  K  f11
mainprogram;  %这是运算出 M  K  f11的系数矩阵的程式
xdot(1)=x(2);
xdot(2)=inv(M)*f11-inv(M)*K*x(1);   
xdot=[xdot(1),xdot(2)];

再来主程式的部分
global M K f11

[t,x]=ode45['inoran',[0,5],[0,0]]

显示错误码   我也知道我这样写会有问题  ,不过我又说不出个所以然
我真的是很初学,身边没啥人可以问,想请问各位高手
我的问题出在那边,亦或者说这种形式的ODE我该怎么撰写？
我试了好久,领悟力真的很差,先谢谢各位！！！！！
小弟一鞠躬！！！！！真的很想把他学好

[ 本帖最后由 ChaChing 于 2010-3-21 23:44 编辑 ]

inoran

各位大侠高手好:
小弟是初学者,我有一微分方程式[M]y''+[k]y=[F]sin(t)
[M]  ,[K]:6*6矩阵       [F]:6*1矩阵　　
我看了书上有   　　二维一阶微分的范例,　与一维二阶微分的范例
现在我卡到我不会用二维（或多维）二阶微分的写法

写进去都会说我矩阵维度有问题,　我是先写个程式计算出M K 等等的系数(以下我先把sin(t)当成1)

然后先将方程式降阶  function  dy=test(t,x)
                                                 global M K F
                                   dx=[x(2);-inv(M)*K*x(1)+inv(M)*F]      %觉得问题在这边,我是参考一阶二维的写法,

                                  主程式：
                                 clear;
                                 global M K F
                                 tinterv=[0 5];
                                 yinit=[2 1];
                                 [t1,y1]=ode113('js',tinterv,yinit);
==========================
　　　　　　　　　　以上我若把参数矩阵改成1*1是一定可以跑出答案
                                  我想请问类似这种系统的ＯＤＥ,也就二阶微分,系数是用矩阵的话,该怎么写
                                  我查了过去的文章,有人提到过,不过他好像是已经会解,只是在找更快速的方法
                                  谢谢各位！！！！！  小弟初学也是第一次来这边,如果有没注意到的礼节也还请各位见谅
                                 谢谢！！<(_ _)>

[ 本帖最后由 ChaChing 于 2010-3-21 23:49 编辑 ]

ChaChing

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

inoran

感謝樓上！！
我再研究一下這位高手的寫法！
希望我也能用他的觀念  寫出解線性couple的程式來！！
感謝感謝！！！！

inoran

请问:原程式码如下
clear all
n=3; F=[25;24;20]; m1=31.2; m2=31.2; m3=31.2;
k1=67.51; k2=67.51; k3=67.51; c1=0.01; c2=0.01;  c3=0.01;
M=[m1,0,0;0,m2,0;0,0,m3];
B=[c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3];
K=[k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3];
DL=inline('[x(n+1:end,1); inv(M)*(F-B*x(1:n,1)-K*x(1:n,1))]',...
          't','x','flag','n','M','K','F','B');
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[t,x]=ode45(DL,[0,3],rand(n,1),options,n,M,K,F,B);  
plot(t,x(:,1:n))
=============

我改成2*2去试试看,改了系数矩阵
F=[25;24]; M=[m1,0;0,m2];
B=[c1+c2,-c2;-c2,c2+c3];
K=[k1+k2,-k2;-k2,k2+k3];
其他的东西我都没改,可是程式不能run  原程式是OK低
小弟的盲点出在哪  谢谢!!

[ 本帖最后由 ChaChing 于 2010-3-22 16:12 编辑 ]

inoran

另外說說我的感想
我來在台灣,在台灣的論壇或是BBS也發問過
不過幾乎找不到會的人,也許還是有啦,只是沒在上論壇
這邊強的真的很多,我也拜讀了許多文章,收穫良多!!

此外我問的問題我也想到可以解的方法
不過我是把矩陣拆開一條一條代換來算
速度上當然很慢也很冗長

謝謝這裡給我的新思維與學習!也謝謝這裡的板友不厭其煩的回答
真的大感謝!!

daivder

回复 1 # inoran 的帖子

请问你的问题解决了么？我要解的方程和你的一样，正愁呢，哪位高手能给我解决一下