ODE solver- can't define a variable used in equation based on previous time steps output

I have a second order differential equation set up using ODE45 to solve for displacement and velocity, but I want to put in a set of statements that define one of the variables in the equation depending on the values from the previous time step:

%The if statements to define the variable I think should be like this...
if (z(2) < 0)
   P_int(round(t/0.1),1)= -5;
elseif   (z(2) > 0)
   P_int(round(t/0.1),1)= 5;
    else
   P_int(round(t/0.1),1)= 0;
end;

%and the general set up of the differential equations
z_dot(1)= z(2);     %velocity of membrane
z_dot(2)= 1/mass*((...+...- P_int +..);
z_dot=z_dot';

But when I check the output of P_int it only gets saves the last entry.

Any help would be much appreciated.

[Star Strider]

What do you want ‘P_int’ to be or do? Is it an array or a function?

What do you want to do with z_dot(2)? The line in your Question:

z_dot(2)= 1/mass*((...+...- P_int +..);

doesn’t convey any useful information.

It’s best to not ‘sketch out’ code. Details are important.

[Jeff]

Sorry I tried to chop it down to keep it neat but made it unclear. P_int is a parameter I want to define depending on the value of z(2), so it would be an array of values for each time step. When z(2)>0 I want to define P_int1 for a certain value, and a different value when z(2)<0 and also for z(2)=0.

 here are my two scripts:

%%%%%%%%%%%%

clear all;
global F1_interp F2_interp F3_interp F4_interp V01 dV1 m01 m1 p01 T0 rho;

t0= 0.1;
tf= 84.5;
tspan= t0:0.1:tf;
%z0= [0;0; 0;0; 0;0; 0;0];      %Initial values
z0= [0;0];

options= odeset('reltol',1e-6,'abstol',1e-6);
[t,z]= ode45(@membraneMotionInterp, tspan, z0);

p01= 113207 + 101325;   %hydrostatic plus atmospheric
rho= 1.2;

subplot(3,1,1);                     %Forces
%plot(t,F1_interp(:,1),'k')
plot(t,F1_interp(:,1),'k')
xlabel('Time (s)');
ylabel('Force (N)');
title('Acting Force');

subplot(3,1,2);                     %Displacement
%plot(t,z(:,1),'k')
plot(t,z(:,1),'k')
xlabel('Time (s)');
ylabel('Translation (m)');
title('Position of Membrane');

subplot(3,1,3);                     %Velocity
%plot(t,z(:,2),'k')
plot(t,z(:,2),'k')
xlabel('Time (s)');
ylabel('Velocity (m/s)');
title('Velocity of Membrane');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function z_dot= membraneMotionInterp(t,z)
%ode23 uses current time and displacement z as input. motion.m calculates
%and returns the rates of change of z, or the velocity of for those
%conditions used by assignment_1

global F1 F2 F3 F4 n m t_irreg F1_irreg F2_irreg F3_irreg F4_irreg t_reg F1_interp F2_interp F3_interp F4_interp  F_P_mean m_WEC n_valves A_valve D_valve rho K_L A T0 V01 p01 dV1 m01 m1 p1 V0_H m0_H V0_L m0_L V0_T m0_T;

F1= csvread('IrregularForcesTest.csv', 1, 0, [1, 0, 1240, 5]);    %reads 6 columns, starting from second row

t_irreg= F1(:,1);
F1_irreg= F1(:,6);
t_reg= 0.1:0.1:84.5;

F1_interp=interp1(t_irreg,F1_irreg,t_reg,'linear')';

A= 16.96;                   %SA of membrane (m^2)
m_WEC= 500000;              %system mass (not a calculated value)
K_L= 30;                    %minor losses coefficient
rho= 1.2;                   %air density
D_valve= 0.2;               %diameter of valve
A_valve= pi*D_valve^2/4;    %area of open valve
n_valves= 4;               %no. valves
F_P_mean= F1_interp(1,1);       %mean water pressure force

P_L= -1000;
P_H= 1000;
if (z(2) < 0)
   P_int(round(t/0.1),1)= -5;
elseif   (z(2) > 0)
   P_int(round(t/0.1),1)= 5;
    else
   P_int(round(t/0.1),1)= 0;
end;

z_dot(1)= z(2);     %velocity of membrane
z_dot(2)= 1/m_WEC*((F1_interp(round(t/0.1),1)-F_P_mean-P_int1-50000*z(1))-n_valves*A_valve*0.5*K_L*rho*(z(2)/A_valve).^2);
z_dot = z_dot';
%%%%%%%%%%%%%%%%%%%%%

from what I can see I need to make an event function instead of having the if statements so the integration stops when z(2) hits zero, but I'm having trouble understanding how to incorporate the variable P_int1 values in this so I was curious if there is another method. Thanks

[Jan]

It is still not clear, why you define P_int as a vector and use round(t/0.1) as index.

There is no reason to apply the time consuming import of the CSV file in each iteration, when the contents is store in the global variable F1.

[Jeff]

I'm fairly new to Matlab so that's just me not understanding properly, from what I can tell I need to use round(t/0.1) in the state space equations to read the correct F1 index for each time step as defined in tspan, I assumed I would also need to define this when creating P_int as an array.

Didn't even think to use CSVread before the loop starts, thankyou!

Answer 1

[Jan]

It is not meaningful to request the value "of the last time step". The function to be integrated is evaluated several times for each time step and as well the step-size control as the event detection can reject steps and reduce the step width. Therefore "previous time step" is not well defined.

Btw. with the shown code your function to be integrated is not smooth. This is a conflict with the specifications of Matlabs ODE45, see http://www.mathworks.com/matlabcentral/answers/59582#answer_72047 .

Answer 2

[Star Strider]

I’m not certain how ‘P_int’ interacts with time, but you can considerably simplify it, completely eliminating the if block:

P_int = @(z2) 5*sign(z2);
q = P_int([-1 0 1])         % Test ‘P_int’

You can simplify it further by eliminating the function expression of it, but the anonymous function form has utility you may want.

[Jeff]

Thanks, but when I check the values of P_int it's still only showing as one value instead of switching depending on z(2). Does the use of ODE solver allow for parameters to change like this?

[Star Strider]

My pleasure!

The ODE solvers allow parameters to change (but it prefers them to change so as not to introduce abrupt discontinuities). If you also run my ‘q’ line, you will see that it works correctly. (I’d not have posted it otherwise.)

What still mystifies me is the relationship ‘P_int’ has with time. It’s not a function of time, and I have no idea what you’re doing by defining it as such. To the best of my knowledge, in the ODE solvers, time is a discrete, single scalar value at each step. It is not a vector.

Also, as Jan Simon notes, if ‘P_int’ changes your equations in a stepwise fashion during your integration, it introduces discontinuities that make your ODE function not smooth. This causes the solver to behave strangely in some situations (see the link).

Answer 3

[Jan]

There is a substantial misconception in your code. Defining P_int as a vector is free of sense. The values of tspan do not matter inside the function to be integrated and there are much more function evaluations than only to the times of tspan. Using round(t/0.1) as an index looks strange and I cannot believe that this has a connection to the mathematical definition of the problem.

The initial time t0 of your problem is 0.1 . Then round(t/0.1) is 0, such that you cannot use it as an index of the vector P_int at all. For a correct code, t0 must be greater than 5. But when the correctness of the code depends on the input values, a programming error is very likely. Notice that if your code works, P_int would be a vector of zeros except for the last element. But in this line:

z_dot(2)= 1/m_WEC*((F1_interp(round(t/0.1),1)-F_P_mean-P_int1 ...

the left-hand-side is a scalar and the right-hand-side is a vector - except for the fact that "P_int1" is not defined, but I assume this is a typo and "P_int" is meant.

The options are set by odeset(), but you did not provide them to the integrator.

So finally I think, that this code cannot be fixed, but should be rewritten from scratch. It is impossible to guess what you want to achieve by P_int, when we see not working code only. Therefore I recommend to update the question by adding the mathematical formula you want to implement.

[Jeff]

OK thanks, looks like I've misunderstood quite a bit. I'll look into how to use an array of values as a forcing function for my differential equations and then try again.