matlab - 使用runge kutta的數值方法的2次二次ODE繪圖

Question

給出以下等式。

我知道我需要打破2二階常微分方程的爲4個一階ODE的。這是我所擁有的。我首先介紹了新的變量u，並且在圖片的底部我寫了我使用ODE45的matlab函數。

現在的問題是，即時通訊應該得到一個拋物線形圖（藍線），但是這不是我所得到。

我已經通過我的代碼沒有結果去了一千倍。可以檢測我的功能中的任何錯誤？

主程序

global g H R alfa 
alfa=pi/2; 
g = 20.0; 
R = 1; 
H=2.3; 
k = 0:0.01:2; 
[T,Y] = ode45(@fspace,k,[H 0 0 0]); 
plot(T,Y(:,1)) 
hold on 


fi= 0:2*pi/60:2*pi; 
xx =R*cos(fi); 
yy =R*sin(fi); 
plot(xx,yy) 

函數f

function f = fspace(x,u) 
global g R H alfa G 
G=(g*R.^2)./((R+H).^2); 
f = [u(2) G*cos(alfa)-g*((R.^2)/u(1).^2)+u(1)*u(4)^2 u(4) (G*sin(alfa)-2*u(2)*u(4))/u(1)]; 

Answer 1

我覺得你的問題在於這些線路：

fi= 0:2*pi/60:2*pi; 
xx =R*cos(fi); 
yy =R*sin(fi); 
plot(xx,yy) 

您正在密謀半徑的圓210。 phi從ode求解器解決方案的一部分，所以你應該改爲：

plot(R*cos(Y(:,3)),R*sin(Y(:,3))) 

但總會給你半徑R的圈子，從來沒有一個拋物線。或者拋物線是指plot(T,Y(:,1))。

方程和代碼看起來是正確的，據我所知。用Y(:,3)替代phi的定義給出了基本相同的圖，但分辨率較低。正如我所說，通過繪製yy與xx，你總會得到一個圓圈。你需要澄清什麼是拋物線應該引用。