2013-10-17 29 views
0

我的代碼給了我零到處爲我的解決方案向量,但我不知道爲什麼。我已經將一個耦合的二階ODE分解爲四階一階ODE。MATLAB ODE解算器給0的地方

我有我的函數定義爲xp.m

function zprime = f(t,z) 
a = 1; 
b = 1; 
c = 1.5; 

zprime = zeros(4,1); 

zprime(1) = z(2); 
zprime(2) = -a*z(1) + b*(z(3) - z(1)); 
zprime(3) = z(4); 
zprime(4) = -c*(z(3) - z(1)); 
end 

我運行它在使用下面的命令MATLAB:

[t,z] = ode45('xp',[1,100],[0 0 0 0]); 

,因爲我的初始條件都爲0。難道我最初的條件給0解決方案還是別的什麼?當我改變集成電路時,解決方案會改變,正如預期的那樣。

感謝

回答

2

您的實際問題與Matlab或編程有關的數學更多。如果你插入0到你的函數f,你會立即看到它除了零之外不會有其他答案。您應該查看equilibrium pointsfixed points。即使一個平衡不穩定(想象一下山頂),一個完全沒有干擾的狀態(外部輸入,數值錯誤)將永遠保持平衡。如果要使用微分方程,那麼知道如何找到平衡點以及如何計算在這些點上評估的雅可比矩陣以確定系統屬性是很好的。如果您在這方面還有其他問題,我建議您去math.stackexchange.com

另外,您正在使用舊方案來調用您的集成函數。你也可以傳入你的參數。定義你的集成功能,在您的主要功能還是作爲一個單獨的函數文件的一個子功能(相同的文件名函數名–你會想比f其他東西)

function zprime = f(t,z,a,b,c) 
zprime(1,1) = z(2); 
zprime(2,1) = -a*z(1) + b*(z(3) - z(1)); 
zprime(3,1) = z(4); 
zprime(4,1) = -c*(z(3) - z(1)); 

然後,在你的主要功能,致電

a = 1; 
b = 1; 
c = 1.5; 
[t,z] = ode45(@(t,z)f(t,z,a,b,c),[1 100],[0 0 0 0]); 
4

爲您的特定情況下,I.C.s z_0 = [0,0,0,0],溶液應平穩,與z_out = [0,0,0,0]值。只需看看你的功能,當你插入z = z_0並運行你的ODE解算器。

zprime(1) = z(2);      % ---> 0 
zprime(2) = -a*z(1) + b*(z(3) - z(1)); % ---> -a*(0) + b*(0) 
zprime(3) = z(4);      % ---> 0 
zprime(4) = -c*(z(3) - z(1));   % ---> -c*(0-0) = 0 

請記住數值解的一般前提。你需要一個初始條件,並將其輸入到導數的公式中。這告訴你正在解決的函數的斜率。您可以使用該斜率來確定將來某個時間(或附近位置或類似位置)的函數值,然後將其反饋回微分公式並重新開始。

不同方法(前進/後退歐拉,RK4,...)之間唯一的主要區別是用於確定當前位置的斜率的方法。