Using a matrix of solution elements outside of ode45

Question

Duncan McIntosh 2016-12-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/316357-using-a-matrix-of-solution-elements-outside-of-ode45

编辑： Walter Roberson 2016-12-12

I am working on a problem for a falling parachute where the forces are proportional to velocity squared. I am using ode45 which calls a function which includes the differential equations of motion. ode45 outputs to the main program time, distance, and velocity. Within the function, I have built a matrix which includes time, distance, velocity, glide path angle, and acceleration using persistent variables. I would like to use this complete matrix after ode45 has finished in the main program but I only get the variables for the last time period and not the whole matrix. How do I make the whole matrix from the function inside the ode45 loop available outside ode45 in the main program?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mischa Kim 2016-12-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/316357-using-a-matrix-of-solution-elements-outside-of-ode45#answer_246749

Duncan, how about simply re-building the matrix after the ode45-call? I assume all those variables are dependent on time, distance, and velocity.

For a more detailed analysis, please attach your code.

2 个评论
显示无隐藏无

Duncan McIntosh 2016-12-11

编辑：Walter Roberson 2016-12-12

在 MATLAB Online 中打开

function [T,T1,Y,Y1]=solution_module;
data_module;
global ACCEL r_o_enu;
v_0=0;a_0=0; alt=r_o_enu(3); % Initial Conditions for unsteady state ode. Parachute being deployed.
[T,Y]=ode45('usm_module',[0:1:5],[v_0;a_0]);
V=Y(6,2); A=ACCEL;  % Initial Conditions (velocity, acceleration) for steady state ode.  Parachute deployed.
[T1,Y1]=ode45('ssm_module',[5:1:450],[V;A])
end
$$$$$$$$$$$$$$$$$$$
function [F]=usm_module(t,x);
data_module;
global W_pl S_bar_pl S_bar_c b e r_o_enu g ACCEL;
alt=r_o_enu(3)+x(1); [rho] = rho_module(alt);
AR=(b.^2)./S_bar_c;  m=W_pl./g;
A=.1.*rho.*S_bar_c.*t;
B=(rho./2).*(S_bar_c.*((.01+(.19./5).*t)+(.2.*t).^2./(pi.*e.*AR))+S_bar_pl.*2);
gama=-atan(B./A).*(180./pi);
k=sqrt(A.^2+B.^2)./(m.*((sind(gama)).^2));
F1=[x(2)];
F2=[k*(x(2)^2)-g];
F=[F1;F2];
ACCEL=[F2];
end
$$$$$$$$$$$$$$$$$$$$$
function [FF]=ssm_module(t,x);
data_module;
global W_pl S_bar_pl S_bar_c b e r_o_enu g;
F1=zeros(500,1);  F2=zeros(500,1); 
alt=r_o_enu(3)+x(1); [rho] = rho_module(alt);
AR=(b.^2)./S_bar_c; m=W_pl./g; t_d=5;
A=.1.*rho.*S_bar_c.*t_d;
B=(rho./2).*(S_bar_c.*((.01+(.19./5).*t_d)+(.2.*t_d).^2./(pi.*e.*AR))+S_bar_pl.*2);
gama=-atan(B./A).*(180./pi);
k=sqrt(A.^2+B.^2)./(m.*((sind(gama)).^2));
FF1=[x(2)];
FF2=[k*(x(2)^2)-g];
FF=[FF1;FF2];
end
$$$$$$$$$$$$$$$$
function data_module;  % The Data Module contains all of the fixed parameters for all the other
modules.  Sets sets these parameters to Global Variables.
global r W_pl S_bar_pl Cd_pl Cs_pl S_bar_c b Cd_c_o Cl_c_o Cd_c_i Cl_c_i e dt_infl lat_o long_o r_o_enu v_ac_atm g;
r=20925670; % (f), radius of the earth
W_pl=100; % (#f), weight of the payload
S_bar_pl=3; % (f^2), characteristic area of payload
Cd_pl=2; % Payload Drag Coefficient 
Cs_pl=.5; % Payload Side Force Coefficient
S_bar_c=20; % (f^2), characteristic area of canopy
b=10; % (f), canopy wing span
Cd_c_o=.01; % Canopy Drag Coefficient, Initial
Cl_c_o=0; % Canopy Lift Coefficinet, Initial
Cd_c_i=.2; % Canopy Drag Coefficient, Inflated
Cl_c_i=1; % Canopy Lift Coefficient, Inflated
e=1; % Oswald Efficiency Factor
dt_infl=5; % (s), Inflation Time
lat_o=44.8831; % (degrees),Initial Latitude (N)
long_o=-68.6719; % (degrees), Initial Longitude (W)
r_o_enu=[0 0 15000]; % (f), Initial position vector (from aircraft)
v_ac_atm=[200 100 0]; % (f/s), Aircraft Initial velocity
g=32.174; % (f/s^2), earth gravitation constant
end
$$$$$$$$$$$$$$$$$$$$$$$$$$$$
function [rho] = rho_module(alt);  % (slug/f^3) Converts geometric to geopotential altitude and then calculates density. Returns density and geometric altitude.
global r g;
temp_sl=518.69;  % (degR), sea level temperature 
rho_sl=.002377; % (slug/f^3), sea level density
a=-.00356; % (degR/f), standard temperature lapse rate
R=1716.49; % (f #f/slug degR)
% Converts geometric to geopotential altitude
alt_gp = (r./(r+alt)).*alt; % (ft), geopotential altitude
% Calculates density using geopotential altitude
temp=temp_sl+a.*alt_gp; % Calculates standard temperature
rho=rho_sl.*((temp./temp_sl).^(-((g./(a.*R))+1))); % Calculates standard density
end

Duncan McIntosh 2016-12-11

May not be the best approach but I am looking at writing the matrix to Excel from within usm_module and ssm_module and then reading the Excel files from solution_module.

请先登录，再进行评论。