如何調試這個Matlab代碼來模擬滑翔機下降？

Question

我有這個Matlab項目，但由於某種原因，我不能停止考慮它，因爲我無法得到它的工作。

目的：

這是一個MATLAB腳本，將計算的壓力，溫度和正在從萬英尺下降滑翔機的密度的變化。當它下降時，我們希望使用這些新值計算，然後插入一個函數，該函數有4個方程，需要使用ode45以及P T和Rho的新值對每個點進行區分。

下面是主要代碼：

% HouseKeeping: 

clc 
clear all 
close all 

% Constants: 
S = 232;           % ft^2 
Cd0 = 0.02; 
K = 0.07; 
W = 11000;          % lbf 
Cl_max = sqrt(Cd0/K); 
Cd_max = 2*K*Cl_max^2; 
Rho_10000 = .001756;        % slugs/ ft^3 

%Initial conditions: 
t = 0;           % Sec 
x = 0;           % ft 
h = 10000;          % ft 
v = sqrt((2*W)/(Rho_10000*S*Cl_max));   %ft/s 
gamma = -Cd_max/Cl_max; 

% Find Endurance: 
V_RD= sqrt((2*W)/(S* Rho_10000* sqrt(3*Cd0/K))); 
RD= V_RD/((sqrt(3*Cd0/K))/(2*Cd0)) ;   % ft/s 
Endurcance= h/RD;        % 958.3515 sec 

% Sea Level values: 
TSL = 518.69;         % Rankine 
PSL = 2116.199414;        % lbs/ft^2 
RhoSL = 0.0023769;       % slugs/ft^3 

while h > 0 

    tspan = [t t+1]; 
    i=1; 
    X = [x;h;v;gamma;Rho_10000]; 
    Time(i)= t; 

    % Calling ODE45: 
    [F] = ode45(@ D,tspan, X) 

    % Hight Varying Parameters: 
    T = TSL - 0.00356616*h; 
    P = (1.137193514E-11)*(T)^5.2560613; 
    Rho = (RhoSL * TSL/PASL)*(P/T); 
    a = 49.0214 * (T)^.5; 

    H_Del(i) = (-Cd_max/Cl_max)*(plotted_x(i))+10000; 
    V_Del(i) = sqrt((2*W)/(Rho*S*Cl_max)); 
    Gamma_Del(i) = -Cd_max/Cl_max; 
    i= i+1; 

    X = [ x ; H_Del(i); V_Del(i); Gamma_Del(i); Rho]; 

end 



% Plots: 

%1 
figure (1) 
plot(F(:,1),'-r',F(:,3),'-b') 
title('Velocity vs Distance'); 
xlabel('x (ft)'); 
ylabel('v (ft/s)'); 


%2 
Figure (2) 
plot(F(:,1),'-r',F(:,2),'-b') 
title('Altitude vs Distance '); 
xlabel('x (ft)'); 
ylabel('h (ft)'); 


%3 
figure (3) 
plot(F(:,1),'-r',F(:,4),'-b') 
title('Gamma vs Distance'); 
xlabel('x (ft)'); 
ylabel('Gamma (rad)'); 

%4 
figure (4) 
plot(t,'-r',F(:,3),'-b') 
title('Velocity vs Time'); 
xlabel(' t (s)'); 
ylabel('v (ft/s)'); 

%5 
figure (5) 
plot(t,'-r',F(:,1),'-b') 
title(' Distance vs Time '); 
xlabel('t (s)'); 
ylabel('x (ft)'); 


%6 
figure (6) 
plot(t,'-r',F(:,4),'-b') 
title('Gamma vs Time '); 
xlabel('t (s)'); 
ylabel('Gamma (rad)'); 

%7 
figure (7) 
plot(t,'-r',F(:,3),'-b') 
title('Velocity vs Time'); 
xlabel('t (s)'); 
ylabel('V (ft/s)'); 

這裏是功能

function [F] = D(X) 

% Constants: 
S = 232;           % ft^2 
Cd0 = 0.02; 
K = 0.07; 
W = 11000;          % lbf 
Cl_max = sqrt(Cd0/K); 
Cd_max = 2*K*Cl_max^2; 
Rho_10000 = .001756;        % slugs/ ft^3 

% Initial conditions: 
t = 0;           % Sec 
x = 0;           % ft 
h = 10000;          % ft 
v = sqrt((2*W)/(Rho_10000*S*Cl_max));   % ft/s 
gamma = -Cd_max/Cl_max; 
g= 32.2;           % ft/s^2 

% Equations: 
X_Pr = X(3)*cos(X(4)); 
H_Pr = X(3)*sin(X(4)); 
V_Pr = (-32.2./W)*(0.5*X(5)*Cd_max*S*X(3)^2 + W*sin(X(4))); 
G_Pr = (32.2./(W*X(3)))*(0.5*X(5)*Cl_max*S*X(3)^2 - W*cos(X(4))); 

F = [X_Pr;H_Pr;V_Pr;G_Pr]; 

我不是很好用Matlab的，但我做了我最好的！我去找我的教授尋求幫助，但他們說他們太忙了。我甚至跟我認識的每一位資深人士說話，他們都說他們不知道該怎麼做。我的下一個項目將盡快分配，我認爲如果我不能做到這一點，那麼我就無法做下一個項目。

Answer 1

你的代碼產生以下錯誤：

Error using main>D

Too many input arguments.

這意味着ODE45嘗試了太多的輸入參數來調用您提供的功能D。您應該檢查所需odefun格式ode45文檔中：dydt = odefun(t,y)

所以，你應該的D你的函數聲明改爲

function [F] = D(t, X) 

這應該解決您的第一個問題，但下面的錯誤彈出：

D returns a vector of length 4, but the length of initial conditions vector is 5. The vector returned by D and the initial conditions vector must have the same number of elements.

同樣，您應該檢查ode45文檔。你的函數應該返回所有輸入變量的衍生物X：F = dX/dt。您還應該返回第五個元素Rho_10000的派生值。

接下來，我得到了一些關於未定義變量的錯誤，如PASL。可能是因爲你沒有發佈完整的代碼。

除了錯誤，你應該再次檢查你的代碼。你寫了一個無限循環while h > 0。循環中永遠不會更改h，也不會在循環中使用ode45的輸出。此外，您總是在循環開始時覆蓋您的i和X值，這可能不是您想要的值。

這不是您的問題的完整答案，但我希望您能夠繼續併發布較小的，明確定義的問題，而不是一個很難完全回答的大問題。