Nested Function is not working correctly with ode45

5 views (last 30 days)
I am trying to use a nested function which uses multiple if conditions to represent time periods and then I integrated using ode45 but it results in a incorrect graph I have no idea why. I am aiming to integrate between 4 time periods one between 0 and tb1, tb1 to tb+tc, tc+tb1 to tb2+tc+tb1 and greater than tb2+tc+t
This is my code:
function EGA215_1234567_A1_Q2
m01=600;
mf1=350;
T1=6500;
Isp1=270;
m02=100;
n2=1.6;
me2=0.5;
Isp2=290;
g0=9.81;
T2=Isp2*g0*me2;
tb1=101.87;
tc=10;
tb2=75;
%boost stage 1
c1=Isp1*g0;
me1=T1/c1;
tspan1=[0 300];
v0=0;
[t,v]=ode45(@rocket_eqn,tspan1,v0);
figure(1)
plot(t,v)
function vdot = rocket_eqn(t, v)
if t<=tb1
vdot=6500/(m01-me1*t)-g0;
elseif (tb1<=t)<=tb1+tc
vdot=-g0;
elseif (tb1+tc<=t)<=tb1+tc+tb2
vdot=T2/(m02-me2*t)-g0;
elseif tb2+tc+tb1<=t
vdot=-g0;
end
end
end

Answers (2)

Walter Roberson
Walter Roberson on 26 Mar 2023
In MATLAB, (A<=B)<=C does not mean to test whether B is between A and C. Instead, it means to test whether A<=B and if so emit 1 and otherwise emit 0, and then to compare that 0 or 1 to C.
EGA215_1234567_A1_Q2
function EGA215_1234567_A1_Q2
m01=600;
mf1=350;
T1=6500;
Isp1=270;
m02=100;
n2=1.6;
me2=0.5;
Isp2=290;
g0=9.81;
T2=Isp2*g0*me2;
tb1=101.87;
tc=10;
tb2=75;
%boost stage 1
c1=Isp1*g0;
me1=T1/c1;
tspan1=[0 300];
v0=0;
[t,v]=ode45(@rocket_eqn,tspan1,v0);
figure(1)
plot(t,v)
function vdot = rocket_eqn(t, v)
if t<=tb1
vdot=6500/(m01-me1*t)-g0;
elseif t<=tb1+tc
vdot=-g0;
elseif t<=tb1+tc+tb2
vdot=T2/(m02-me2*t)-g0;
elseif tb2+tc+tb1<=t
vdot=-g0;
else
vdot = nan;
end
end
end
  1 Comment
Walter Roberson
Walter Roberson on 26 Mar 2023
Edited: Walter Roberson on 26 Mar 2023
However, when you use ode45, the mathematics of Runge-Kutta requires that the first two derivatives of your function are continuous. In the great majority of cases, when people use if inside an ode function, they have not carefully ensured that the first and second derivatives are continuous at the boundaries.
If you are lucky, when you use if inside of an ode function, you will receive an error message indicating that ode45 was unable to integrate over a singularity.
If you are not lucky, then instead you will simply get the wrong result and not realize that it is the wrong result.
The above plot is not the right result.
You should study the ballode example to see how to use event functions.

Sign in to comment.


Torsten
Torsten on 26 Mar 2023
Edited: Torsten on 26 Mar 2023
You can also get an analytical solution for v since your piecewise equation is easily integrated.
EGA215_1234567_A1_Q2
function EGA215_1234567_A1_Q2
m01=600;
mf1=350;
T1=6500;
Isp1=270;
m02=100;
n2=1.6;
me2=0.5;
Isp2=290;
g0=9.81;
T2=Isp2*g0*me2;
tb1=101.87;
tc=10;
tb2=75;
%boost stage 1
c1=Isp1*g0;
me1=T1/c1;
%tspan1=[0 300];
v0=0;
tspan1 = [0 tb1];
iflag = 1;
[t1,v1]=ode45(@(t,y)rocket_eqn(t,y,iflag),tspan1,v0);
tspan2 = [tb1 tb1+tc];
iflag = 2;
[t2,v2]=ode45(@(t,y)rocket_eqn(t,y,iflag),tspan2,v1(end));
tspan3 = [tb1+tc tb1+tc+tb2];
iflag = 3;
[t3,v3]=ode45(@(t,y)rocket_eqn(t,y,iflag),tspan3,v2(end));
tspan4 = [tb1+tc+tb2 300];
iflag = 4;
[t4,v4]=ode45(@(t,y)rocket_eqn(t,y,iflag),tspan4,v3(end));
t = [t1;t2;t3;t4];
v = [v1;v2;v3;v4];
figure(1)
plot(t,v)
function vdot = rocket_eqn(t, v,iflag)
if iflag==1
vdot=6500/(m01-me1*t)-g0;
elseif iflag==2
vdot=-g0;
elseif iflag==3
vdot=T2/(m02-me2*t)-g0;
elseif iflag==4
vdot=-g0;
end
end
end

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!