[問題]請問一振動問題

看板Physics作者 (貓與鋼琴)時間15年前 (2010/04/18 17:58), 編輯推噓5(5推 0噓 79→)
留言84則, 8人參與, 6年前最新討論串1/2 (看更多)
大家好： 我有一懸臂樑,自由端受一簡諧力,簡化後用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)
→
04/18 19:56, , 1F
我覺得起始條件y''(0)不等於零
04/18 19:56, 1F
→
04/18 21:10, , 2F
我自己是想說初始加速度＝0
04/18 21:10, 2F
→
04/18 21:10, , 3F
不過我帶入ODE 找出的y''(0)也是0
04/18 21:10, 3F
→
04/18 21:10, , 4F
其實我也是一直覺得我的初始條件有問題
04/18 21:10, 4F
→
04/18 21:11, , 5F
y(0)=y'(0)應該是沒錯才對 y''(0)=0我就不太確定
04/18 21:11, 5F
推
04/18 21:43, , 6F
這樣想好了
04/18 21:43, 6F
→
04/18 21:43, , 7F
妳竟然有外力了
04/18 21:43, 7F
→
04/18 21:43, , 8F
代表有加速度 y''(0) 了
04/18 21:43, 8F
→
04/18 21:43, , 9F
故不為零才對
04/18 21:43, 9F
→
04/18 21:50, , 10F
可是我的外力是sin(10t) t=0時 外力項=0
04/18 21:50, 10F
先謝謝回答！ 回代原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)
→
04/18 22:47, , 11F
所以算出來只有y1不一樣嗎 還是y1~y6兩個方法算都不一樣
04/18 22:47, 11F
→
04/18 23:04, , 12F
都不一樣
04/18 23:04, 12F
→
04/18 23:05, , 13F
因為y2~y6對系統影響很小 所以可以不用太在意
04/18 23:05, 13F
→
04/18 23:05, , 14F
我還是覺得問題在初始條件,但是y''(0)=0這樣做出來
04/18 23:05, 14F
→
04/18 23:06, , 15F
也沒錯 所以我找不到原因＠＠
04/18 23:06, 15F
→
04/18 23:07, , 16F
FFT轉成頻域圖的話 取log大概差了兩個order
04/18 23:07, 16F
→
04/18 23:07, , 17F
看時域圖的話 理論上是要在x=0上下震盪
04/18 23:07, 17F
→
04/18 23:07, , 18F
往上偏移的話那就怪了
04/18 23:07, 18F
→
04/18 23:12, , 19F
哪表示物體不是通過水平線的往復振動
04/18 23:12, 19F
推
04/18 23:32, , 20F
我覺得問題可能在於 原來的ODE是有物理意義的 也就是牛頓定律
04/18 23:32, 20F
→
04/18 23:33, , 21F
而你對整個ODE微分一次 在數學上是成立的 但是微分之後整條式
04/18 23:33, 21F
→
04/18 23:34, , 22F
子代表的物理意義卻是?? 這可能是出問題的地方
04/18 23:34, 22F
謝謝E大回答,這也是我之前思考過的地方, 我會再認真思考看看！ 謝謝！ ※ 編輯: inoran54787 來自: 122.116.13.191 (04/18 23:39)
→
04/19 07:41, , 23F
程式寫錯了? 你有把解代回去檢查嗎?
04/19 07:41, 23F
→
04/19 07:42, , 24F
你的M11 k11 f1 和兩個方程式的解 y1 分別是甚麼?
04/19 07:42, 24F
程式我是這樣寫的(三階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)
→
04/19 15:01, , 25F
看起來沒錯. 怪怪的..我隨便給個M f k, 跑出來一樣
04/19 15:01, 25F
→
04/19 15:04, , 26F
只不過..這微分方程應該不需要跑程式. 手解就可以
04/19 15:04, 26F
→
04/19 20:50, , 27F
@@
04/19 20:50, 27F
→
04/19 20:51, , 28F
就是因為是18*18 然後還有耦合項
04/19 20:51, 28F
→
04/19 20:52, , 29F
看不懂..但是我兩個方法跑出來結果相同.
04/19 20:52, 29F
→
04/19 20:52, , 30F
這邊我只是先把耦合項拿掉
04/19 20:52, 30F
→
04/19 20:52, , 31F
手解可以吧..
04/19 20:52, 31F
→
04/19 20:52, , 32F
C大有興趣的話我可以給你我的程式
04/19 20:52, 32F
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)
→
04/19 20:53, , 33F
y = A sin(wt) + Bsin(10t). A B 都是 vector
04/19 20:53, 33F
推
04/19 20:58, , 34F
基本上沒辦法用手解拉 那麼大的維度手解會死人
04/19 20:58, 34F
→
04/19 21:00, , 35F
不可以嗎? 我對ODE是不太熟練. 但是 應該是上述解吧.
04/19 21:00, 35F
→
04/19 21:01, , 36F
之後矩陣運算如果不想手解, 交給mathematica 或
04/19 21:01, 36F
→
04/19 21:01, , 37F
matlab 都可以.
04/19 21:01, 37F
→
04/19 21:01, , 38F
h學長來了＠＠
04/19 21:01, 38F
三階主程式 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)
→
04/19 21:02, , 39F
if t<=0.1 這邊有點怪..why? 這變造成不同吧
04/19 21:02, 39F
→
04/19 21:06, , 40F
改成小於也可啦 我是想在0～0.1有個SIN的脈衝波
04/19 21:06, 40F
→
04/19 21:07, , 41F
不是啊..你給了 t<=0.1. 這樣在這範圍微分不是cos.
04/19 21:07, 41F
→
04/19 21:07, , 42F
這邊就會讓你真正的啟始值 y(0.1) 不同
04/19 21:07, 42F
→
04/19 21:08, , 43F
t<=0.1不就是 0~0.1的意思？
04/19 21:08, 43F
推
04/19 21:09, , 44F
你的M矩陣裏面 怎麼會有-0? 那是實0還是虛0?
04/19 21:09, 44F
→
04/19 21:09, , 45F
所以我的初始直還是在y(t=0)
04/19 21:09, 45F
→
04/19 21:10, , 46F
都是時數 只是在寫個程式算出來的東西
04/19 21:10, 46F
→
04/19 21:10, , 47F
我直接複製貼上去低
04/19 21:10, 47F
→
04/19 21:10, , 48F
不影響 兩者比較過後的值,因為兩者帶的參數都一樣
04/19 21:10, 48F
→
04/19 21:11, , 49F
-0 是負數趨近於０
04/19 21:11, 49F
推
04/19 21:12, , 50F
你手動把他改成0 再跑看看吧
04/19 21:12, 50F
→
04/19 21:12, , 51F
好 我試試看＠＠
04/19 21:12, 51F
→
04/19 21:13, , 52F
我不是這意思..你 三階ODE 的 cos 怎麼寫?
04/19 21:13, 52F
※ 編輯: inoran54787 來自: 140.115.66.117 (04/19 21:13)
→
04/19 21:13, , 53F
在 t<=0.1 的部分你怎麼寫?
04/19 21:13, 53F
→
04/19 21:14, , 54F
你是有改過內容嗎? 不然就是我剛剛看錯了 :p
04/19 21:14, 54F
→
04/19 21:16, , 55F
沒有也
04/19 21:16, 55F
→
04/19 21:17, , 56F
我只是改f11 變成f 我原程式寫的是f11
04/19 21:17, 56F
→
04/19 21:18, , 57F
oh..那是我看錯了
04/19 21:18, 57F
→
04/19 21:18, , 58F
三階的我寫在上面了 二階的也有＠＠ 勞駕了！
04/19 21:18, 58F
→
04/19 21:19, , 59F
哈哈 沒關係 程式常常會讓人看到眼花
04/19 21:19, 59F
→
04/19 21:38, , 60F
所以..你的 y 是在 t<0.1開始就不同了?
04/19 21:38, 60F
→
04/19 21:38, , 61F
或是 t>0.1 才開始不同?
04/19 21:38, 61F
→
04/19 21:40, , 62F
不 程式跑出來的y是t=0~5
04/19 21:40, 62F
→
04/19 21:40, , 63F
我是0～5秒出來的解整個都不同
04/19 21:40, 63F
→
04/19 21:40, , 64F
我知道啊..我是想知道你的 y 從何時開始不同.
04/19 21:40, 64F
→
04/19 21:40, , 65F
程式計算是把t回代ODE計算出y
04/19 21:40, 65F
→
04/19 21:41, , 66F
0秒開始就不一樣
04/19 21:41, 66F
→
04/19 21:41, , 67F
所以連 t<0.1 也不同嗎?
04/19 21:41, 67F
→
04/19 21:41, , 68F
oh...
04/19 21:41, 68F
→
04/19 21:41, , 69F
震盪的位移變的很小
04/19 21:41, 69F
→
04/19 21:41, , 70F
YAP
04/19 21:41, 70F
→
04/19 21:42, , 71F
二階做出來在-1.5~1.5震盪（10^-4）
04/19 21:42, 71F
→
04/19 21:42, , 72F
如果是0.1之後開始不同, 就稍微能理解.
04/19 21:42, 72F
→
04/19 21:42, , 73F
三階做出來在3~3.5震盪（10^-4）
04/19 21:42, 73F
→
04/19 21:43, , 74F
這個方法不行我也只能放棄T.T
04/19 21:43, 74F
→
04/19 21:50, , 75F
我先回去重新確定我的IC
04/19 21:50, 75F
→
04/19 21:50, , 76F
因為我把0,-0都當成0來計算 所以我的y''(0)=0
04/19 21:50, 76F
→
04/19 21:51, , 77F
這是我不夠謹慎的地方,C大也謝謝你
04/19 21:51, 77F
→
04/19 21:52, , 78F
待我確認完 我明天在回復 感謝鼎力相助!
04/19 21:52, 78F
其實我主要的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)
→
08/13 15:34, , 79F
我覺得起始條件y''( https://noxiv.com
08/13 15:34, 79F
→
09/17 10:14, , 80F
待我確認完 我明天在回 https://daxiv.com
09/17 10:14, 80F
→
09/17 13:30, , 81F
y = A sin(w https://daxiv.com
09/17 13:30, 81F
→
11/09 11:09, , 82F
程式寫錯了? 你有把解 https://muxiv.com
11/09 11:09, 82F
→
01/02 14:23, 7年前 , 83F
震盪的位移變的很小 https://muxiv.com
01/02 14:23, 83F
→
07/06 21:58, 6年前 , 84F
oh..那是我看錯了 http://yofuk.com
07/06 21:58, 84F
‣ 更多 inoran54787 的文章
文章代碼(AID): #1BojVIrT (Physics)
‣ 更多 inoran54787 的文章
文章代碼(AID): #1BojVIrT (Physics)