声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2088|回复: 6

[编程技巧] 高阶系数矩阵的二阶微分方程

[复制链接]
发表于 2010-3-20 23:39 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
各位高手好 :
小弟是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 编辑 ]
回复
分享到:

使用道具 举报

 楼主| 发表于 2010-3-21 22:40 | 显示全部楼层

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

各位大侠高手好:
小弟是初学者,我有一微分方程式[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是一定可以跑出答案
                                  我想请问类似这种系统的ODE,也就二阶微分,系数是用矩阵的话,该怎么写
                                  我查了过去的文章,有人提到过,不过他好像是已经会解,只是在找更快速的方法
                                  谢谢各位!!!!!  小弟初学也是第一次来这边,如果有没注意到的礼节也还请各位见谅
                                 谢谢!!<(_ _)>

[ 本帖最后由 ChaChing 于 2010-3-21 23:49 编辑 ]
发表于 2010-3-22 00:06 | 显示全部楼层
 楼主| 发表于 2010-3-22 00:30 | 显示全部楼层
感謝樓上!!
我再研究一下這位高手的寫法!
希望我也能用他的觀念  寫出解線性couple的程式來!!
感謝感謝!!!!
 楼主| 发表于 2010-3-22 02:12 | 显示全部楼层
请问:原程式码如下
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 编辑 ]
 楼主| 发表于 2010-3-22 02:20 | 显示全部楼层
另外說說我的感想
我來在台灣,在台灣的論壇或是BBS也發問過
不過幾乎找不到會的人,也許還是有啦,只是沒在上論壇
這邊強的真的很多,我也拜讀了許多文章,收穫良多!!

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

謝謝這裡給我的新思維與學習!也謝謝這裡的板友不厭其煩的回答
真的大感謝!!
发表于 2011-5-31 15:33 | 显示全部楼层
回复 1 # inoran 的帖子

请问你的问题解决了么?我要解的方程和你的一样,正愁呢,哪位高手能给我解决一下
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 09:49 , Processed in 0.057219 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表