Solving nonlinear second order diferential equation
3 views (last 30 days)
Show older comments
Dados iniciais
clear
Massa_foguete_seco = 3000; % Kg
Massa_combustivel_inicial = 3000; % kg
Empuxo = 300000; % Newton (sentido de y)
Consumo_combustivel = 100; % Kg/segundo
C_arrasto = 50; % Newton/(metros/segundo)
Gravidade = 9.81; % Metros/segundo^2 (sentido contrario de y)
rho = 2; % Kg/metros cubicos
Area_secao = 15; % Metros quadrados
Peso do Foguete.
syms y(t) t
Massa_combustivel = Massa_combustivel_inicial - Consumo_combustivel * t;
Massa_foguete = Massa_combustivel + Massa_foguete_seco;
Peso_foguete = Massa_foguete * -Gravidade;
Cálculo da Aceleração
Forca_arrasto = -0.5 * C_arrasto * rho * Area_secao * diff(y,t)^2;
Forca_empuxo = Empuxo;
Peso = Peso_foguete;
Forca_resultante = Forca_arrasto + Forca_empuxo + Peso;
Equações do movimento
Aceleracao = Forca_resultante/ Massa_foguete
ode = diff(y,t,2) == Aceleracao;
dy = diff(y);
cond1 = y(0) == 0;
cond2 = dy(0) == 0;
cond = [cond1 cond2];
Altitude = dsolve(ode,cond)
Velocidade = diff(Altitude);
Aceleracao = diff(Altitude,2);
0 Comments
Answers (1)
Stephan
on 15 Oct 2020
%% Dados iniciais
clear
Massa_foguete_seco = 3000; % Kg
Massa_combustivel_inicial = 3000; % kg
Empuxo = 300000; % Newton (sentido de y)
Consumo_combustivel = 100; % Kg/segundo
C_arrasto = 50; % Newton/(metros/segundo)
Gravidade = 9.81; % Metros/segundo^2 (sentido contrario de y)
rho = 2; % Kg/metros cubicos
Area_secao = 15; % Metros quadrados
%% Peso do Foguete.
syms y(t)
Massa_combustivel = Massa_combustivel_inicial - Consumo_combustivel * t;
Massa_foguete = Massa_combustivel + Massa_foguete_seco;
Peso_foguete = Massa_foguete * -Gravidade;
%% Cálculo da Aceleração
Forca_arrasto = -0.5 * C_arrasto * rho * Area_secao * diff(y,t)^2;
Forca_empuxo = Empuxo;
Peso = Peso_foguete;
Forca_resultante = Forca_arrasto + Forca_empuxo + Peso;
%% Equações do movimento
Aceleracao = Forca_resultante/ Massa_foguete;
ode = diff(y,t,2) == Aceleracao
%% Solve numerically
[eqs, vars] = odeToVectorField(ode);
odefun = matlabFunction(eqs,'vars',{'t','Y'})
tspan = [0 5];
init = [0 0];
[t,y] = ode45(odefun,tspan,init);
plot(t,y)
2 Comments
Stephan
on 18 Oct 2020
Edited: Stephan
on 19 Oct 2020
This is a numerical solution. One line represents the velocity, the other the way. Since this is a second order system, which is transformed into a first order system of 2 ode, you get the solution to both, the way and the velocity.
See here also:
To calculate the acceleration from the results use:
%% Dados iniciais
clear
Massa_foguete_seco = 3000; % Kg
Massa_combustivel_inicial = 3000; % kg
Empuxo = 300000; % Newton (sentido de y)
Consumo_combustivel = 100; % Kg/segundo
C_arrasto = 50; % Newton/(metros/segundo)
Gravidade = 9.81; % Metros/segundo^2 (sentido contrario de y)
rho = 2; % Kg/metros cubicos
Area_secao = 15; % Metros quadrados
%% Peso do Foguete.
syms y(t)
Massa_combustivel = Massa_combustivel_inicial - Consumo_combustivel * t;
Massa_foguete = Massa_combustivel + Massa_foguete_seco;
Peso_foguete = Massa_foguete * -Gravidade;
%% Cálculo da Aceleração
Forca_arrasto = -0.5 * C_arrasto * rho * Area_secao * diff(y,t)^2;
Forca_empuxo = Empuxo;
Peso = Peso_foguete;
Forca_resultante = Forca_arrasto + Forca_empuxo + Peso;
%% Equações do movimento
Aceleracao_num = Forca_resultante/ Massa_foguete;
ode = diff(y,t,2) == Aceleracao_num;
%% Solve numerically
[eqs, vars] = odeToVectorField(ode);
odefun = matlabFunction(eqs,'vars',{'t','Y'});
tspan = [0 5];
init = [0 0];
[t,y] = ode45(odefun,tspan,init);
% Calculate acceleration from numeric solution - Create ab function handle
% that represents the acceleration ode, but now we know the numeric
% solution for the ode system, so we can use the velocity informations
% needed in the original ode.
Aceleracao_num = @(t,Y2) -(981*t - 750*Y2.^2 + 241140)./(100*t - 6000);
% Calculate those values based on the results from ode45 and also save them
% in y-vector, column 3
y(:,3) = Aceleracao_num(t,y(:,2));
% plot the whole system, with y, y_dot and y_doubledot
plot(t,y)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!