2012-11-21 22 views
1

我正在使用一個簡單的if循環來更改我的ode腳本中的參數值。這是我寫的一個示例腳本,它展示了同樣的問題。因此,首先其作品的版本:循環中的時間相關參數值在頌歌求解器中僅適用於某些t的值

function aah = al(t,x) 

if (t>10000 && t<10300) 
    ab = [0; 150]; 
else 
    ab = [150; 0]; 
end 

aah = [ab]; 

這可以使用

t = [0:1:10400]; 
x0 = [0,0]; 
[t,x] = ode23tb(@al, t,x0); 

運行,並與

plot(t,x(:,1)) 
plot(t,x(:,2)) 

確定這就是好版可視化。現在如果你所做的只是改變t到

t = [0:1:12000]; 

整件事情都爆炸了。你可能會認爲它只是matlab的平均出圖,但它不是,因爲如果你看一下

​​

答案應該是在這兩種情況下相同,因爲代碼並沒有改變。但是這第二個版本輸出0,這是錯誤的!

究竟是怎麼回事?任何人有想法?

非常感謝你的幫助

回答

3

你的函數是恆定的(除了10000<噸< 10300),因此內部求解器開始通過解決系統具有非常大的時間步長的總時間的10%默認。 (在自適應ODE求解器中,如果系統是常數,高階和低階解決方案將給出相同的解,並且(估計)誤差將爲零,因此求解器假定當前時間步長足夠好)。看看你是否只給了tspan兩個元素,開始和結束時間。

t = [0 12000]; 

通常tspan不會影響解算器的內部時間步長。求解器用它們的內部時間步長求解系統,然後插入用戶給出的tspan。因此,如果內部時間步驟不幸「跳過」時間間隔[10000,10300],解算器將不知道間隔。

所以你最好設置的最大步長,相對小於300

options = odeset('MaxStep', 10); 
[t, x] = ode23tb(@al, t, x0, options); 

如果你不想用小的步長全部的時間來解決(如果你「知道」當函數是不是恆定的),你應該分開解決。

t1 = [0 9990]; 
t2 = [9990 10310]; 
t3 = [10310 12000]; 

[T1, x1] = ode23tb(@al, t1, x0); 
[T2, x2] = ode23tb(@al, t2, x1(end,:)); 
[T3, x3] = ode23tb(@al, t3, x2(end,:)); 

T = [T1; T2(2:end); T3(2:end)]; 
x = [x1; x2(2:end,:); x3(2:end,:)]; 
+0

實際上,您可以將步長設置爲300.大多數頌歌算法將「及時回溯」以糾正其可能錯過的快速更改。只要一步落入敏感區域,你就會好起來 – Rasman

+0

絕對美妙。非常有意義,非常感謝你。 – user1792403

相關問題