Re: [問題] 非線性ODE求解已回收
經過前輩指正，
寫完程式碼如下:
function CaseII
global a b c m e
a=0.4;
b=15055*10^(-6);
c=523;
m=120;
e=0.8;
[t x] = ode15s(@CaseII20100306,[0 500],[0 298])
plot(t,x)
function f = CaseII20100306(t,x)
global a b c m e
syms z g
y =(4422-(((80-10*cos(t)-(70^2-(10*sin(t))^2)^0.5)*201)));
z =(4422-(((80-10*cos(g)-(70^2-(10*sin(g))^2)^0.5)*201)));
f = a.*(x./y).*subs(diff(z),t)+b.*(x-c)*((x+m)/(y.*x^1.5))^e*(y./4+1105.87);
==============================================================================
但結果為發散，無法收歛，
於是去除熱傳項，
f = a.*(x./y).*subs(diff(z),t) %+b.*(x-c)*((x+m)/(y.*x^1.5))^e*(y./4+1105.87);
計算結果仍無法呈現週期函數。
直接將 y 微分得 z 帶入運算式，
程式碼改為:
====================================================
function CaseII
global a b c m e
a=0.4;
b=15055*10^(-6);
c=523;
m=120;
e=0.8;
[t x] = ode15s(@CaseII20100306,[0 500],[0 298])
plot(t,x)
function f = CaseII20100306(t,x)
global a b c m e
syms z g
y =(4422-(((80-10*cos(t)-(70^2-(10*sin(t))^2)^0.5)*201)));
z=-2010*sin(t)-20100/(4900-100*sin(t)^2)^(1/2)*sin(t)*cos(t);
f = -a*x*z/y;
=============================================================================
結果仍為發散...
與手算的結果不符合，
應該要是週期函數才對，
請問為什麼會這樣呢?
我有特別畫出 z/y 圖形為週期函數，
照理說 plot(t,x)
應該也會是週期函數才對，
但是算出來數值卻會越來越大，
請問是哪個環節有出錯嗎?
感謝!!!
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 122.118.2.88
※ 編輯: YangJimmy 來自: 122.118.2.88 (03/09 01:09)
推
05/05 01:28, , 1F
05/05 01:28, 1F
討論串 (同標題文章)
完整討論串 (本文為第 3 之 4 篇)：
問題
1
3