Borrar filtros
Borrar filtros

Solve differential equation Systems with ODE45

1 visualización (últimos 30 días)
Brendan Görres
Brendan Görres el 22 de Oct. de 2018
Comentada: Brendan Görres el 23 de Oct. de 2018
I want to solve a differental equation system with ODE45.
close all;clc;
%basic planet/moon Parametres%
re=6378e3;
g0=9.81;
rho0=1.225;
sh=7500;
%parametres related to the vehicle%
m0=68e3;
Cdrag=0.5;
diameter=5;
area= pi*((diameter/2)^2);
Adrag= area*Cdrag;
Isp=390;
thrust= 933910;
R=15;
%initial calculations%
ceff= Isp*g0;
propflow= thrust/ceff;
mf=m0/R;
mp= m0-mf;
tburn=mp/propflow;
%Functions%
syms h;
rhoh=rho0*exp(-h/sh);
gh=g0*((re/re+h)^2);
r(h)=re+h;
%Pitchover altitude mass and velocity%
syms x
y0=89.85;
hpo=-(g0/2)*(((m0-x)/propflow)^2)+(ceff/propflow)*(x*log(x/m0)+m0-x)==130;
mpo=vpasolve(hpo,x);
vpo=-g0*((m0-mpo)/propflow)-ceff*log(mpo/m0);
tpo=(m0-mpo)/propflow;
%inital conditions%
tRange=[0 tburn];
Y0=[y0;vpo;hpo;0;mpo];
%Solve ODE§
[tSol,YSol]=ode45(@AscentProblem,tRange,Y0);
%Define Functions%
function dYdt= AscentProblem(t,Y)
y=Y(1);
v=Y(2);
h=Y(3);
x=Y(4);
m=Y(5);
dydt=((v(t)/(re+h(t)))-(g0/v(t)).*((re/re+h(t)).^2)).*cos(y(t));
dvdt=(propflow*ceff/m(t))-(Adrag* rhoh.*(v(t).^2))/2*m(t)- gh .*sin(y(t));
dhdt=v(t).*sin(y(t));
dxdt=(re/(re+h(t))).*v(t).*cos(y(t));
dmdt=-propflow;
%Define dYdt%
dYdt=[dydt;dvdt;dhdt;dxdt;dmdt];
end
Array indices must be positive integers or logical values.
Error in sym/subsref (line 870) R_tilde = builtin('subsref',L_tilde,Idx);
Error in APode45>AscentProblem (line 50) dydt=((v(t)/(re+h(t)))-(g0/v(t)).*((re/re+h(t)).^2)).*cos(y(t));
Error in odearguments (line 90) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in APode45 (line 40) [tSol,YSol]=ode45(@AscentProblem,tRange,Y0);
PLease help, I cannot find my mistake.

Respuesta aceptada

Torsten
Torsten el 22 de Oct. de 2018
Y0=[y0;vpo;hpo;0;mpo];
must be a vector of doubles.
Check the type of every entry of this vector and modify your code accordingly.
  3 comentarios
Torsten
Torsten el 22 de Oct. de 2018
Replace v(t), h(t) y(t) and m(t) by v,h,y and m.
Further make all the constants you use in AscentProblem accessible there (g0, re, propflow,...)
Brendan Görres
Brendan Görres el 22 de Oct. de 2018
Now MatLab gives me these errors.
Error in odearguments (line 90) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in APode45 (line 43) [tSol,YSol]=ode45(@AscentProblem,tRange,Y0);

Iniciar sesión para comentar.

Más respuestas (1)

Steven Lord
Steven Lord el 22 de Oct. de 2018
The variables you define inside your ODE function as each containing one element of the Y input with which ode45 calls your ODE function (y, v, h, x, and m) are scalars. You shouldn't be trying to index into them using the variable t. t is unlikely to be the exact value 1, meaning it is unlikely to be a valid index into a scalar. Just use y, v, h, x, and m.

Categorías

Más información sobre Programming en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by