ode45 command matlab help

1000*dV/dt+cd*1.225*V^2*A+200-T=0 this is the ode i have to currently solve for assignment due on the 11th of this month where T is defined by the for loop added below. Heres the question "Write a function into which you pass the end time and time step size that returns the time and speed array for a car initially at rest, speed equals 0.0 m/s, is accelerated along a runway " cd=0.5 and A=2

for t=0:time_step:time_end
if t<1
    T=10000*t
elseif 1<=t<=50
    T=10000
end
end

[Steven Lord]

This command doesn't do what you think it does.

1<=t<=50

That is equivalent to:

(1<=t) <=50

Since (1 <= t) is either true or false, which when compared to a double precision number are treated as 1 or 0, this expression is always true no matter what t is.

(1<=t) & (t<=50)

This will return true if t is greater than or equal to 1 and less than or equal to 50.

With that out of the way, show us what you've tried to do to solve this problem and ask a specific question about where you're having difficulty and we may be able to offer some guidance.

[Abdullaziz Errayes]

i tried to create a separate function with the ODE in it then called the function using the ode command but theres one problem with that the T is undefined. i dont really know how to implement the ode function in one script file that has the for loop with the if statement in order to define T

function dvdt=ode1(t,v)
A=2;
dvdt=(T-200-(cd*1.225*v^2*A))/1000;
end

and this is my initial script

for t=0:time_step:time_end
if t<1
    T=10000*t
elseif (1<=t) & (t<=50)
    T=10000
end
end 
tspan=0:time_step:time_end;
[t,v]=ode45(@ode1 ,tspan,0)

plot(t,v)

[Torsten]

function main
time_step = ...;
time_end = ...;
tspan=0:time_step:time_end;
[t,v]=ode45(@ode1 ,tspan,0);  
plot(t,v)

function dvdt=ode1(t,v)
A=2;
T=10000*(t*(t<1) + (t>=1)*(t<=50));
cd = ...;
dvdt=(T-200-(cd*1.225*v^2*A))/1000;
end

Best wishes

Torsten.

[Abdullaziz Errayes]

i appreciate you taking the time to help me but if i run the code posted above theres an error message that says "The expression to the left of the equals sign is not a valid target for an assignment." for the time_end=...; and

[Torsten]

Of course, you will have to fill the "..." with numbers :-)

[Abdullaziz Errayes]

yeah i figured i had to as soon as i commented haha now my next problem is even if i do fill it with numbers it doesn't show an error but it doesn't do anything when i run it, am i missing something obvious ahha

[Abdullaziz Errayes]

time_end=input('Input end time');
time_step=input('Input time step size');
function main
time_step =input('Input end time');
time_end= input('Input time step size');
tspan=0:time_step:time_end;
[t,v]=ode45(@ode2 ,tspan,0);  
plot(t,v)
function dvdt=ode2(t,v)
A=2;
T=10000*(t*(t<1) + (t>=1)*(t<=50));
cd = 0.5;
dvdt=(T-200-(cd*1.225*v^2*A))/1000;
end
end

here is my current code but when i run it theres no error but it doesnt do anything

[Jan]

@Abdullaziz Errayes: If you have saved this code to a script file, nothing happens, because the function main() is not called at all. Remove the first two lines and the code will run.

"a function into which you pass the end time and time step size" sounds like your function should take 2 input arguments, and not two input commands.

Answer 1

[Jan]

Remember that Matlab's ODE integrators (as many other tools also) can handle smooth functions only: The automatic step size control gets deeply confused if the function to be integrated is not differentiable. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. If you are lucky, the integrator stops with an error, but without luck, you get a result, which is dominated by the accumulated rounding error, because the step size is reduced until the discontinuity is hidden by the truncation error.

The numerically correct method is to stop the integration at the discontinuity and restart it. If the time is known in advance, e.g. as Tstop:

[t1, y1] = ode45(fcn1, [0, Tstop], y0);
[t2, y2] = ode45(fcn2, [Tstop, Tend], y1(end, :));
t = [t1; t2];
y = [y1; y2];

If the discontinuity depends of a certain position, use an event function to stop the integration and a while loop until t(end) is the wanted final time.

[EDITED] The question of your homework is not precise. Maybe a simple Euler method is sufficient already instead of a smart ODE integrator or symbolic integrations:

function [t, v] = YourFcn(tEnd, dt)
t  = 0:dt:tEnd;
v0 = 0;
cd = 0.5;
A  = 2;

v    = zeros(size(t));
v(1) = v0;   % Actually not needed because it is 0 already
for it = 2:numel(t)
  % if t < 1, T = 10000*t; else, T = 10000; end  -- easier:
  T     = 10000 * min(t(it), 1);
  accel = (T - 200 - (cd * 1.225 * v(it-1) ^ 2 * A)) / 1000;
  v(it) = v(it - 1) + accel * dt;
end
end

[John BG]

Jan,

this is John BG <mailto:[e-mail] [e-mail]>

you went for ode45

then I posted my answer, no need for ode45.

and now you add that 'maybe' (!?) a simple euler method suffices ..

either you need ode45, or you don't, what's going to be?

[Jan]

@John BG:

"you went for ode45" - wrong. The OP Abdullaziz Errayes has mentioned ODE45 and I've explained the common pitfall of integrating a discontinuous function with a step size controlled method.

"either you need ode45, or you don't" - wrong. The OP did not mentioned, if the homework question specified the kind of integration. Then "maybe" a simple Euler method is enough. By the way, the Euler method and the ODE45 approach posted by the OP compute equivalent velocities (in spite of the numerical instability):

The velocity calculated by your suggested symbolic integration looks different. A short test: The final velocity is reached, when the acceleration is 0:

T  = 10000;
cd = 0.5;
A  = 2;
accel = (T - 200 - (cd * 1.225 * v_end ^ 2 * A)) / 1000 = 0  ==>
v_end = sqrt((T - 200) / (cd * 1.225 * A)) = 89.4427

The final v with the Euler method is 89.4393, with ODE45 it is 89.4113, your value is about 810 according to the diagram. I suggest to cross-check your suggestion instead of posting a pointless comment.

Answer 2

[John BG]

Hi Abdullaziz

this is John BG <mailto:[e-mail] [e-mail]>

1.

the ODE is simple and it can be split into 2 parts:

1000*dV/dt+cd*1.225*V^2*A+200    -    T

a.

a basic 1st degree ODE that doesn't really need any solver to solve it, its just the integration of a polynomial:

1000*dV/dt+cd*1.225*V^2*A+200

b.

The Torque applied

- T

2.

How does the Torque and its integral look like?

clear all;close all;clc

time_end=120   % 2 minute
time_step=.1
A=2;
cd = 0.5;

% time_step =input('Input end time');
% time_end= input('Input time step size');

% tspan=[0:time_step:time_end];
t=[0:time_step:time_end];
T=10000*(t.*(t<1) + (t>=1).*(t<=50));
T1=cumsum(T)

yyaxis left
plot(t,T);grid on
ylabel('Torque')
axis([-10 130 -100 10500])

yyaxis right
plot(t,T1)
xlabel('t(sec)')
ylabel('int(Torque)')

title(' Torque and int(Torque)')
axis([-5 125 -100 5e6])

.

.

An example of real Torque is available in page 26 of Balaji Kamalammannan's Master Degree thesis

'Modelling and Simulation of Vehicle Kinematics and Dynamics (WITH MATLAB)'

.

link:

https://hh.diva-portal.org/smash/get/diva2:1115795/FULLTEXT02.pdf

3.

Just split the needed integration

   int(f1 + f2)= inf(f1) + int(f2)

syms dvdt v
fv=1/1e3*(cd*1.225*v^2+200)
v_without_T=int(fv)
v_without_T =
(49*v^3)/240000 + v/5

Then

v=v_without_T + T1

T1 being cumsum(T)

4.

Since the integral of the Torque is already null at v=0, let's find the integration constant to meet the initial condition

v(0)=0

v=(49*v^2)/80000 + 1/5+T1  + kv0

kv0=-.5-T1(0)

5.

Solving

kv1=-.5    % to set v(0)=0
fv3=@(u) (49*u^2)/80000 + 1/5+T1k+ kv1 - u

v=0
for k=1:1:numel(t)
    T1k=T1(k);
    fv3=@(u) (49*u^2)/80000 + 1/5+T1k+ kv1 - u;
    d2=eval(solve(fv3,0));
    v=[v real(d2(1))];
end

6.

plotting

figure;plot(v);grid on
axis([0 100 -100 900])

.

.

As expected the dominant term of the velocity, with the equation you use, is the Torque.

Torque is roughly proportional to fuel (or battery charge) consumption.

But your equation is for an ideal model, with the vehicle on flat terrain. In reality, only with certain minimum downhill angle, and not even with brand new bearings and everything working well, air drag, friction continuously reduces gained velocity thus further torque is required to keep speed constant.

Abdulaziz

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance for time and attention

John BG

<mailto:[e-mail] [e-mail]>