If statements in ode45 function

i am using the following signal as an input in a matlab code that i have made.I have put it inside the ode45 function file but the results that i am getting don't make sense. I have tried using the fixed step ode4 solver and the results are what i expect to be, but i don't want to use a fixed step solver as i don't trust them to produce accurate results.

t = linspace(0,100,10000);
Fmax = 10000;
ts1 = 1;
ts2 = 2;
ts3= 0;
ts4=2;
tt = ts1+ts2+ts3+ts4;
n = max(size(t));
Fn = zeros(1,n);
for i = 1:n
    B = floor(t(i)/(tt));
    if t(i) <= (B)*tt+ts1
        Fn(i) =Fmax.*((t(i)-B*tt)/ts1); 
    elseif t(i) > (B)*tt+ts1 && t(i) <= (B)*tt+ts1+ts2
        Fn(i) = Fmax;
    elseif t(i) > (B)*tt+ts1+ts2 && t(i) <= (B)*tt+ts1+ts2+ts3
        Fn(i) = -Fmax.*(t(i)-(ts2+ts1)-B*tt)/abs(ts3-ts2)+Fmax;
    elseif t(i) > (B)*tt+ts1+ts2+ts3 && t(i) < (B)*tt+ts1+ts2+ts3+ts4
        Fn(i) = 0;
    end
end
plot(t,Fn)

Of course inside ode function i have taken out the ''for'' part and the i.

Thanks

Answer 1

[Star Strider]

The MATLAB ODE solvers are adaptive-step, not fixed-step solvers. Everything I’ve read about them indicates they’re the most robust algorithms available.

I don’t understand what you’re doing with your code. As a general rule, if you have different conditions (variable values, and such) for different times in your integration, the best strategy is to stop the integration, use the last values for the initial conditions of the next integration, then restart the integration. Repeat this for multiple discontinuities. The ODE solvers do not integrate well across discontinuities.

[theo13]

What you are saying makes sense to me and explains the results that i am getting.Is there any way to locate these discontinuities?And can you give an example of how to stop the integration and start with the last values as initial conditions?Thanks

[Star Strider]

My pleasure.

As an example of the way I would approach this, my first ‘tspan’ vector would be:

tspan = [0 (B)*tt+ts1];

then in your first ODE function, since this is dependent on time, you would define ‘Fn’ as:

Fn = Fmax.*((t-B*tt)/ts1);

then run your first ode function for that time:

[t{1},x{1}] = ode45(@odefun_1, tspan, ic_vct);

the for the second time interval:

tspan = [(B)*tt+ts1,  (B)*tt+ts1+ts2];

with ‘odefun_2’ having:

Fn = Fmax;

and:

[t{2},x{2}] = ode45(@odefun_2, tspan, x{1}(end,:));

and so for the rest.

This is just a sketch, so you will need to experiment with it. I used a cell array for the integrated output to make those data more tractable since they are over different time intervals, and so it would be easier to work with them rather than concatenating them all into one matrix. You can concatenate them easily later if you want, but I would keep them separate at least initially.

There are obviously different ways to approach this. For example, you can have the various values of ‘Fn’ in a switch-case structure in one ODE function, and pass an integer to your ODE function as a separate parameter to switch them for each time interval.

[theo13]

Thanks a lot i will try what you suggested.

[Star Strider]

My pleasure!

[theo13]

I am having some difficulties in the application of the above suggestion.I have pasted my approach below

U0 = 90;
ts1 = 1;
ts2 = 2;
ts3= 0;
ts4=1;
tt_1 = ts1+ts2+ts3+ts4;
t = linspace(0,300,10000);
n = max(size(t));
Fn  = zeros(1,n);

for i = 1:n
    A = floor(t(i)/(tt_1));
    if t(i) <= (A)*tt_1+ts1
        [t_1,x_1]=ode45(@Thermal_function,[0 (A)*tt_1+ts1],initialvalues);
    elseif t(i) > (A)*tt_1+ts1 && t(i) <= (A)*tt_1+ts1+ts2
        [t_2,x_2]=ode45(@Thermal_functiona,[(A)*tt_1+ts1  (A)*tt_1+ts1+ts2],x_1(end,:));
    elseif t(i) > (A)*tt_1+ts1+ts2 && t(i) <= (A)*tt_1+ts1+ts2+ts3
        [t_3,x_3]=ode45(@Thermal_functionb,[(A)*tt_1+ts1+ts2 (A)*tt_1+ts1+ts2+ts3],x_2(end,:));
    elseif t(i) > (A)*tt_1+ts1+ts2+ts3 && t(i) < (A)*tt_1+ts1+ts2+ts3+ts4
       [t_4,x_4]=ode45(@Thermal_functionc,[(A)*tt_1+ts1+ts2+ts3 (A)*tt_1+ts1+ts2+ts3+ts4],x_3(end,:));
    end
end

So the first function file function will be dx=Thermal_function(t_1,x_1) dx=zeros(4,1);% because i am solving 4 differential equations Fn =Fmax.*((t_1-A*tt_1)/ts1); etc. for the other three function files

Does it make any sense ?is what you were suggesting ?

[Star Strider]

That is essentially what I intended, but not exactly.

I don’t know what you are doing with ‘A’, so I didn’t use it. I would just use separate non-overlapping time intervals or vectors instead, something like this:

[t_1,x_1]=ode45(@Thermal_function,[0 (A)*tt_1+ts1],initialvalues);
[t_2,x_2]=ode45(@Thermal_functiona,[(A)*tt_1+ts1  (A)*tt_1+ts1+ts2],x_1(end,:));
[t_3,x_3]=ode45(@Thermal_functionb,[(A)*tt_1+ts1+ts2 (A)*tt_1+ts1+ts2+ts3],x_2(end,:));
[t_4,x_4]=ode45(@Thermal_functionc,[(A)*tt_1+ts1+ts2+ts3 (A)*tt_1+ts1+ts2+ts3+ts4],x_3(end,:));

You don’t need a loop considering how you have created your ODE functions and time intervals. I cannot run your code (and I don’t understand your time intervals), but this modification should work. You have to remove ‘A’ because I removed the loop, since you don’t need it. Just substitute whatever ‘A’ was supposed to calculate in each ode45 call.

Let me know how it goes. I’ll help as much as I can.

[theo13]

A is the same as B in the first post that i made.I will explain what i want to do in order to have a better insight on what i want to achieve. The first code that i have posted its a trapezoidal ''like'' pulse that continuously repeat it self until the end of the simulation.This trapezoidal signal represent a speed that is going to be used as an input inside the ode45 function.A lot of values are going to be calulated based on that speed an at the end four first order differential equations are going to be solved. Correct me if i am wrong, but the above way that you suggested is going to work only for one pulse and not for the whole duration of the simulation. That's way i have put it in a loop, does it make any sense?and again thanks your help

[Star Strider]

My pleasure.

O.K. Now I understand the loop. The code without the loop would just run once. (Running it multiple times will only produce the results from the final iteration, since you’re not saving the results of the previous iterations. If you want to save all of them, there are a few different ways to do it, depending on what you want as the result.)

The ‘A’ may be necessary in that instance, but the if block would still not be, providing I understand correctly what you’re doing.

[theo13]

Is it possible to have the signal without the use of the if?( i need to inform you that in terms of matlab skills i am still a noob!) Also can you saw me an example of how to save the results of the previous iterations?

[Star Strider]

I don’t know what the if block is doing with respect to your signal or simulation. I can’t figure it out by looking at your code.

Don’t worry — we’re all still learning.

The way I would save the results of each iteration would be:

TFc{i,:} = [t_1,x_1;  t_2,x_2;  t_3,x_3;  t_4,x_4];

I tested a version of this with a loop, and it worked. Put it as the last statement in the loop. It vertically concatenates the time vectors and integration results in a matrix, then saves that matrix to a particular cell in ‘TFc’.