Runge-Kutta system equation solver for a system of equations

8 次查看（过去 30 天）
David Wonus
David Wonus 2023-4-16
评论： David Wonus 2023-4-16
Ran in:
Trying to solve a system of equations using the runge-kutta method for estimating the velocity (v) of free fall at various altitudes (z) where z = 0:40000 km
CD = coefficient of dragg
m = mass
A = cross sectional area
The first equation is ds/dt = v the second is dv/dt = gravity - ((CD/m) * ((1/2)*density*(v^2)*A))
where gravity and density are both functions of (z)
this is what I have:
Re = 6371000
Re = 6371000
denEstimated = -0.000000000000037549844448166229302696435310363*z.^3 + 0.0000000035489894922510675713650514186874*z.^2 - 0.00011259491846037372752859645474999*z + 1.2180219634561690877916362296673
Unrecognized function or variable 'z'.
gravity = 9.81 * ((Re.^2)/((z)+Re).^2)
CDFree = 0.735
AFree = 1
m = 120
% setup
dydt = @(t, y) [y(2); gravity-((CDFree/m)*((1/2)*denEstimated(1:401) * (y(2).^2)*AFree))]
y0 = 0:1
t = 0:260
h = 0.5
y = rk4sysn(dydt, t, y0, h)
% Compare results
figure;
hold on;
plot(t, y); plot(t, y, '--', 'linewidth', 2)
function [tp,yp] = rk4sysn(dydt,tspan,y0,h,varargin)
% rk4sys: fourth-order Runge-Kutta for a system of ODEs
% [t,y] = rk4sys(dydt,tspan,y0,h,p1,p2,...): integrates
% a system of ODEs with fourth-order RK method
% input:
% dydt = name of the M-file that evaluates the ODEs
% tspan = [ti, tf]; initial and final times with output
% generated at interval of h, or
% = [t0 t1 ... tf]; specific times where solution output
% y0 = initial values of dependent variables
% h = step size
% p1,p2,... = additional parameters used by dydt
% output:
% tp = vector of independent variable
% yp = vector of solution for dependent variables
if nargin<4,error('at least 4 input arguments required'), end
if any(diff(tspan)<=0),error('tspan not ascending order'), end
n = length(tspan);
ti = tspan(1);tf = tspan(n);
if n == 2
t = (ti:h:tf)'; n = length(t);
if t(n)<tf
t(n+1) = tf;
n = n+1;
end
else
t = tspan;
end
tt = ti; y(1,:) = y0;
np = 1; tp(np) = tt; yp(np,:) = y(1,:);
i=1;
while(1)
tend = t(np+1);
hh = t(np+1) - t(np);
if hh>h,hh = h;end
while(1)
if tt+hh>tend,hh = tend-tt;end
k1 = dydt(tt,y(i,:),varargin{:})'
ymid = y(i,:) + k1.*hh./2;
k2 = dydt(tt+hh/2,ymid,varargin{:})';
ymid = y(i,:) + k2*hh/2;
k3 = dydt(tt+hh/2,ymid,varargin{:})';
yend = y(i,:) + k3*hh;
k4 = dydt(tt+hh,yend,varargin{:})';
phi = (k1+2*(k2+k3)+k4)/6;
y(i+1,:) = y(i,:) + phi*hh;
tt = tt+hh;
i=i+1;
if tt>=tend,break,end
end
np = np+1; tp(np) = tt; yp(np,:) = y(i,:);
if tt>=tf,break,end
end
end

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

回答（1 个）

Walter Roberson
Walter Roberson 2023-4-16
Are you trying to run the ode on a whole series of elevations simultaneously, or does one of your state variables represent the evalation?
gravity = 9.81 * ((Re.^2)/((z)+Re).^2)
If you are trying to run a whole series of elevations simultaneously, then here z would need to be the list of elevations. Your value for Re appears to be Earth radius in metres, so your z would have to be a list of elevations in metres that you wanted to process simultaneously, and you would need to change the / to ./
But if one of your state variables is elevation, then you would need gravity to be a function of an input evalation, not a constant based upon a list of pre-determined elevation. You would need to use an anonymous function, and in the ode you would need to pass the elevation state variable to the anonymous function.
denEstimated = -0.000000000000037549844448166229302696435310363*z.^3 + 0.0000000035489894922510675713650514186874*z.^2 - 0.00011259491846037372752859645474999*z + 1.2180219634561690877916362296673
z is not defined. Again, are you wanting to process a list of pre-defined elevations, or are you wanting to process according to a state variable?
dydt = @(t, y) [y(2); gravity-((CDFree/m)*((1/2)*denEstimated(1:401) * (y(2).^2)*AFree))]
Are you indexing a vector denEstimated at locations 1:401 ? Or are you invoking a function with input [1, 2, 3, ... 401] ? Either way you are probably expecting that the result would be length 401, which would tend to imply that you are running a number of different elevations simultaneously. If denEstimated is made into a function, then are the values 1 to 401 in the proper units for the function? Is that a density estimate for a height in metres like the gravity needs? If so then does it need metres above the surface of the Earth or does it need metres relative to the centre of the Earth ?
  2 个评论
David Wonus
David Wonus 2023-4-16
I have a maxtrix of elevations of length 401 from 0:40000m both gravity and density of air are dependent on the elevations (z) these elevations are relative to the surface such that 0 = sea level
David Wonus
David Wonus 2023-4-16
tspan = [0 260];
y0 = [0; 0]; % initial conditions for s and v
h = 0.5; % step size
[tp,yp] = rk4sys(@myODEFree,tspan,y0,h);
z = yp(:,1)/1000; % extract solution for altitude and convert to km
v = yp(:,2)*3.6; % extract solution for velocity and convert to km/h
plot(z,v)
ylabel('Velocity (km/h)')
xlabel("Altitude Descent (km)")
tspan = [260 560];
y0 = [2567; 0]; % initial conditions for s and v
h = 0.5; % step size
[tp,yp] = rk4sys(@myODEPara,tspan,y0,h);
z = yp(:,1)/1000; % extract solution for altitude and convert to km
v = yp(:,2)*3.6; % extract solution for velocity and convert to km/h
plot(z,v)
ylabel('Velocity (km/h)')
xlabel("Altitude Descent (km)")
function dydt = myODEFree(t,y)
% Constants
CD = 0.735; % Drag coefficient
m = 120; % kg
A = 1; % m^2
z = y(1);
v = y(2);
g = 9.81 * (6371000^2)./((6371000+z).^2);
density = -0.000000000000037549844448166229302696435310363*z.^3 + 0.0000000035489894922510675713650514186874*z.^2 - 0.00011259491846037372752859645474999*z + 1.2180219634561690877916362296673;
dydt = zeros(2,1);
dydt(1) = v;
dydt(2) = g - ((CD/m)*((1/2)*density*v^2*A));
end
% Define a function that evaluates the ODEs
function dydt = myODEPara(t,y)
% Constants
CD = 1.75; % Drag coefficient
m = 120; % kg
A = 25; % m^2
z = y(1);
v = y(2);
g = 9.81 * (6371000^2)./((6371000+z).^2);
density = -0.000000000000037549844448166229302696435310363*z.^3 + 0.0000000035489894922510675713650514186874*z.^2 - 0.00011259491846037372752859645474999*z + 1.2180219634561690877916362296673;
dydt = zeros(2,1);
dydt(1) = v;
dydt(2) = g - ((CD/m)*((1/2)*density*v^2*A));
end
%I am still using the rk4sys code, I no longer have errors but it seems
%like something is wrong as theoretically the velocity should increase
%until reaching terminal velocity until eventually decreasing through the
%atmosphere as air resistance increases with lower altitude

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by