作者inoran54787 (貓與鋼琴)
看板Physics
標題[問題]請問一振動問題
時間Sun Apr 18 17:58:37 2010
大家好：
我有一懸臂樑,自由端受一簡諧力,簡化後用ODE表示如下
初始條件y(0)=y'(0)=0
[M]y''+[k]y=[f]*sin(10*t)
M（6*6） K(6*6) f(6*1)為常數矩陣(因為是連體系統 所以我只取前六個特徵來計算)
我用程式跑出解y1~y6 在物理上表示位移
因為程式需求 所以我對ODE多取一次微分
[M]y'''+[k]y'=[f]*10*cos(10*t)
初始條件y(0)=y'(0)=y''(0)=0
在去跑程式,y1~y6仍為位移
不過前後兩個y1算出來的答案不一樣
也就是原始算出來的位移會在橫軸(y=0)上下震盪
而我多微分一次的解大約在y=2上下震盪
請問這在物理觀點上算正常嗎？,數學上來講的話兩者的解應該是會一樣的才對
小弟好困惑 ,謝謝！！！
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 122.116.13.191
※ 編輯: inoran54787 來自: 122.116.13.191 (04/18 18:00)
1F：→ CHTM:我覺得起始條件y''(0)不等於零 04/18 19:56
2F：→ inoran54787:我自己是想說初始加速度＝0 04/18 21:10
3F：→ inoran54787:不過我帶入ODE 找出的y''(0)也是0 04/18 21:10
4F：→ inoran54787:其實我也是一直覺得我的初始條件有問題 04/18 21:10
5F：→ inoran54787:y(0)=y'(0)應該是沒錯才對 y''(0)=0我就不太確定 04/18 21:11
6F：推 ntust661:這樣想好了 04/18 21:43
7F：→ ntust661:妳竟然有外力了 04/18 21:43
8F：→ ntust661:代表有加速度 y''(0) 了 04/18 21:43
9F：→ ntust661:故不為零才對 04/18 21:43
10F：→ inoran54787:可是我的外力是sin(10t) t=0時 外力項=0 04/18 21:50
先謝謝回答！
回代原ODE y''(0)=-inv[M][k]y(0)+inv[M][f]*sin(0)
y''(0)=0 我自己有點搞混了,找不出原因
※ 編輯: inoran54787 來自: 140.115.66.117 (04/18 21:54)
11F：→ Emcc:所以算出來只有y1不一樣嗎 還是y1~y6兩個方法算都不一樣 04/18 22:47
12F：→ inoran54787:都不一樣 04/18 23:04
13F：→ inoran54787:因為y2~y6對系統影響很小 所以可以不用太在意 04/18 23:05
14F：→ inoran54787:我還是覺得問題在初始條件,但是y''(0)=0這樣做出來 04/18 23:05
15F：→ inoran54787:也沒錯 所以我找不到原因＠＠ 04/18 23:06
16F：→ inoran54787:FFT轉成頻域圖的話 取log大概差了兩個order 04/18 23:07
17F：→ inoran54787:看時域圖的話 理論上是要在x=0上下震盪 04/18 23:07
18F：→ inoran54787:往上偏移的話那就怪了 04/18 23:07
19F：→ inoran54787:哪表示物體不是通過水平線的往復振動 04/18 23:12
20F：推 Emcc:我覺得問題可能在於 原來的ODE是有物理意義的 也就是牛頓定律 04/18 23:32
21F：→ Emcc:而你對整個ODE微分一次 在數學上是成立的 但是微分之後整條式 04/18 23:33
22F：→ Emcc:子代表的物理意義卻是?? 這可能是出問題的地方 04/18 23:34
謝謝E大回答,這也是我之前思考過的地方,
我會再認真思考看看！ 謝謝！
※ 編輯: inoran54787 來自: 122.116.13.191 (04/18 23:39)
23F：→ chungweitw:程式寫錯了? 你有把解代回去檢查嗎? 04/19 07:41
24F：→ chungweitw:你的M11 k11 f1 和兩個方程式的解 y1 分別是甚麼? 04/19 07:42
程式我是這樣寫的(三階ODE)
令x1=y ,x7=y' x13=y'' (不用x1 x2 x3的順序是方便寫程式)
帶入[M]y'''+[k]y'=[f]*10*cos(10*t)
所以y'''=-inv[M][K]y'+inv[M][f]10*cos(10*t)]
降階可寫成
dx1=x7
dx7=x13
dx13=-inv[M][K]x7+inv[M][f]10*cos(10*t)
套成矩陣寫法(MATLAB)
dx =[x(7:12,1)
x(13:18,1)
-inv[M][K]x(7:12,1)+inv[M][f]10*cos(10*t)]
(7:12,1)表示 x7~x12 , 6列1行矩陣
等號左邊dx表示dx1~dx18 18列1行矩陣
以上是18*18矩陣 ,這樣寫基本上我覺得沒問題
M K f寫出來很長 所以我就省略了,因為這裡是物理版,想說我PO程式怪怪的!
※ 編輯: inoran54787 來自: 140.115.66.117 (04/19 12:11)
25F：→ chungweitw:看起來沒錯. 怪怪的..我隨便給個M f k, 跑出來一樣 04/19 15:01
26F：→ chungweitw:只不過..這微分方程應該不需要跑程式. 手解就可以 04/19 15:04
27F：→ inoran54787:@@ 04/19 20:50
28F：→ inoran54787:就是因為是18*18 然後還有耦合項 04/19 20:51
29F：→ chungweitw:看不懂..但是我兩個方法跑出來結果相同. 04/19 20:52
30F：→ inoran54787:這邊我只是先把耦合項拿掉 04/19 20:52
31F：→ chungweitw:手解可以吧.. 04/19 20:52
32F：→ inoran54787:C大有興趣的話我可以給你我的程式 04/19 20:52
main.m
fnction dy=main(t,y)
global M K n f
if t<=0.1
V=sin(10*t)
else
V=0
end
dy=[y(n+1:12,1);-inv(M)*K*y(1:6,1)+inv(M)*f*V]
主程式：
clear all
clc
global M K n f
n=6
M=[0.0324,-0,0,-0,0,-0
0,0.0324,-0,0,-0,0
0,0,0.0324,-0,0,-0
0,0,0,0.0324,-0,0
0,0,0,0,0.0324,-0.0001
-0,-0,0,0,0,0.0324]
K=[185266.167,-141.5699,-3260.0615,-7820.7975,-7835.0081,2155.0833
-141.5699,7263405.7962,1663.5988,3359.6854,2001.6193,-4373.1318
-3260.0615,1663.5988,56978851.3922,77113.5924,73133.927,-31692.6866
-7820.7975,3359.6854,77113.5924,218851300.5359,185151.8109,-51451.9216
-7835.0081,2001.6193,73133.927,185151.8109,597758337.7498,1571.7597
2155.0831,-4373.1317,-31692.6867,-51451.9214,1571.7596,1333597202.0219]
f=[2.7244;-1.9638;2.0016;-1.9999;2.0001;-2.0001]
tinterv=0:0.0002:5
yinit=[0 0 0 0 0 0 0 0 0 0 0 0 ]
options=odeset('RelTol',1e-4)
[t,y]=ode23tb('main',tinterv,yinit,options)
Y=y(:,1)
plot(t,Y);
※ 編輯: inoran54787 來自: 140.115.66.117 (04/19 21:01)
33F：→ chungweitw:y = A sin(wt) + Bsin(10t). A B 都是 vector 04/19 20:53
34F：推 harry901:基本上沒辦法用手解拉 那麼大的維度手解會死人 04/19 20:58
35F：→ chungweitw:不可以嗎? 我對ODE是不太熟練. 但是 應該是上述解吧. 04/19 21:00
36F：→ chungweitw:之後矩陣運算如果不想手解, 交給mathematica 或 04/19 21:01
37F：→ chungweitw:matlab 都可以. 04/19 21:01
38F：→ inoran54787:h學長來了＠＠ 04/19 21:01
三階主程式
fnction dy=main(t,y)
global M K n f
if t<=0.1
V=10*cos(10*t)
else
V=0
end
dy=[y(7:12,1);y(13:18,1);-inv(M)*K*y(7:12,1)+inv(M)*f*V]
ic y(0)=y'(0)=y''(0)=0 如果有興趣的話可以試試看
也謝謝您剛剛的回答 orz
※ 編輯: inoran54787 來自: 140.115.66.117 (04/19 21:05)
39F：→ chungweitw:if t<=0.1 這邊有點怪..why? 這變造成不同吧 04/19 21:02
40F：→ inoran54787:改成小於也可啦 我是想在0～0.1有個SIN的脈衝波 04/19 21:06
41F：→ chungweitw:不是啊..你給了 t<=0.1. 這樣在這範圍微分不是cos. 04/19 21:07
42F：→ chungweitw:這邊就會讓你真正的啟始值 y(0.1) 不同 04/19 21:07
43F：→ inoran54787:t<=0.1不就是 0~0.1的意思？ 04/19 21:08
44F：推 harry901:你的M矩陣裏面 怎麼會有-0? 那是實0還是虛0? 04/19 21:09
45F：→ inoran54787:所以我的初始直還是在y(t=0) 04/19 21:09
46F：→ inoran54787:都是時數 只是在寫個程式算出來的東西 04/19 21:10
47F：→ inoran54787:我直接複製貼上去低 04/19 21:10
48F：→ inoran54787:不影響 兩者比較過後的值,因為兩者帶的參數都一樣 04/19 21:10
49F：→ inoran54787:-0 是負數趨近於０ 04/19 21:11
50F：推 harry901:你手動把他改成0 再跑看看吧 04/19 21:12
51F：→ inoran54787:好 我試試看＠＠ 04/19 21:12
52F：→ chungweitw:我不是這意思..你 三階ODE 的 cos 怎麼寫? 04/19 21:13
※ 編輯: inoran54787 來自: 140.115.66.117 (04/19 21:13)
53F：→ chungweitw:在 t<=0.1 的部分你怎麼寫? 04/19 21:13
54F：→ chungweitw:你是有改過內容嗎? 不然就是我剛剛看錯了 :p 04/19 21:14
55F：→ inoran54787:沒有也 04/19 21:16
56F：→ inoran54787:我只是改f11 變成f 我原程式寫的是f11 04/19 21:17
57F：→ chungweitw:oh..那是我看錯了 04/19 21:18
58F：→ inoran54787:三階的我寫在上面了 二階的也有＠＠ 勞駕了！ 04/19 21:18
59F：→ inoran54787:哈哈 沒關係 程式常常會讓人看到眼花 04/19 21:19
60F：→ chungweitw:所以..你的 y 是在 t<0.1開始就不同了? 04/19 21:38
61F：→ chungweitw:或是 t>0.1 才開始不同? 04/19 21:38
62F：→ inoran54787:不 程式跑出來的y是t=0~5 04/19 21:40
63F：→ inoran54787:我是0～5秒出來的解整個都不同 04/19 21:40
64F：→ chungweitw:我知道啊..我是想知道你的 y 從何時開始不同. 04/19 21:40
65F：→ inoran54787:程式計算是把t回代ODE計算出y 04/19 21:40
66F：→ inoran54787:0秒開始就不一樣 04/19 21:41
67F：→ chungweitw:所以連 t<0.1 也不同嗎? 04/19 21:41
68F：→ chungweitw:oh... 04/19 21:41
69F：→ inoran54787:震盪的位移變的很小 04/19 21:41
70F：→ inoran54787:YAP 04/19 21:41
71F：→ inoran54787:二階做出來在-1.5~1.5震盪（10^-4） 04/19 21:42
72F：→ chungweitw:如果是0.1之後開始不同, 就稍微能理解. 04/19 21:42
73F：→ inoran54787:三階做出來在3~3.5震盪（10^-4） 04/19 21:42
74F：→ inoran54787:這個方法不行我也只能放棄T.T 04/19 21:43
75F：→ inoran54787:我先回去重新確定我的IC 04/19 21:50
76F：→ inoran54787:因為我把0,-0都當成0來計算 所以我的y''(0)=0 04/19 21:50
77F：→ inoran54787:這是我不夠謹慎的地方,C大也謝謝你 04/19 21:51
78F：→ inoran54787:待我確認完 我明天在回復 感謝鼎力相助! 04/19 21:52
其實我主要的ODE是長這樣
[M1]y''+[k1]y=[f1]*sin(10*t)+[c1]∫i(t)dt (1)
[M2]i''+[k2]i'+[c2]i=[c3]y' (2)
解聯立 我現在是想到令∫i(t)dt=x i=x' i'=x'' i''=x'''
假設我已知∫i(0)dt=0
這樣去解 ,(1)的形式就不會去變動到
假使有誤差 也還在容範圍內,本體的自由震動曲線不要跑掉就好
不知這樣可行啦 不過我是打算這樣處理 也謝謝大家的熱心相助
學到不少!
※ 編輯: inoran54787 來自: 140.115.66.117 (04/20 15:21)
※ 編輯: inoran54787 來自: 140.115.67.220 (04/20 20:46)
※ 編輯: inoran54787 來自: 140.115.67.220 (04/20 20:50)