Solve differential equation Systems with ODE45

2 views (last 30 days)
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.

Accepted Answer

Torsten
Torsten on 22 Oct 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 Comments
Torsten
Torsten on 22 Oct 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 on 22 Oct 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);

Sign in to comment.

More Answers (1)

Steven Lord
Steven Lord on 22 Oct 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.

Tags

Community Treasure Hunt

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

Start Hunting!