Problem with solving a system of differential equations using ode45

I am trying to solve a system of differential equations regarding a piston with heat transfer and friction. In this sytem, the friction value is determined by the direction of the motion of the piston. So when the piston is moving to the right, the friction constant is negative, and when it is moving to the left the friction constant is positive. This is fine in the code, however when I solve the equations using ode45 and graph them, the graphs look strange. Below is a velocity vs time graph of the piston.

You can see that extra oscillation later on in the graph. I don't think this oscillation is supposed to happen based on the notes I was given. Another problem I am having is it seems that this oscillation doesn't occur at the same location when I change the time span.

Below are two graphs: one (left) is the graph above but the limits are changed to show the oscillation, the second graph (right) is when the system is only solved between times 116 and 117.

I am trying to analyze whether or not this oscillation is supposed to happen but I don't understand why the graphs change when I change the time span. This also happens when I solve using the time span from 0 to a value closer to 117. Shouldn't the graphs be the same at this location regardless of the time span?

Below is my code for the system of differential equations.

function dydt = odes(~,y)

x = y(1) ;

v = y(2) ;

temp = y(3) ;

temp_p = y(4) ;

%constants

c = .1 ; % t_c / t_h

a = .1 ; %alpha

p_i = 2 ; %initial p

temp_i = 1 ; %initial temp

temp_p_i = 1 ; %initial piston temp

friction_constant = .1 ;

% friction constant based on direction

if y(2) < 0

f = friction_constant ;

else

f = -friction_constant ;

end

p = (temp ./ x) * p_i ;

dxdt = v ;

dvdt = p - 1 + f ;

dtempdt = -((c) * (temp - temp_p)) - ((1 / p_i ) * (2/3) * (p .* v)) ;

dtemp_pdt = (a * c * (temp - temp_p)) + ((a / p_i) * (2/3) * f * abs(v)) ;

dsgdt = ((3/2) * (1 ./ temp) .* dtempdt) + (v ./ x);

dspdt = (1/a) * (3/2) * ((1 ./ temp_p) .* dtemp_pdt) ;

dydt = [dxdt ; dvdt ; dtempdt ; dtemp_pdt ; dsgdt ; dspdt] ;

end

[Jon]

I don't have time at the moment to look at your problem in detail, but one thing you should be careful about is setting the friction force to zero when the velocity is zero. Looking quickly at your code it looks like you may not cover that case in your if statement for computing the friction force

Answer 1

[Star Strider]

I have no idea what the iniitial conditions are suposed to be.

The oscillation might be the result of the if block defining ‘f’, because it creates a non-differentiable step. I usually deal with such transitions by using the tanh function, because it can provide an abrupt transition and still be differentiable. The alternative is to use an event function (see ODE Event Location) to stop the integration each time the piston reverses direction and define a new value for ‘f’.

Here, I implemented that as:

f = -friction_constant*tanh(y(2)*10);

I encourage you to experiment with it to see if it makes a difference in the result.

tspan = [0 150];
ic = zeros(6,1)+0.5;                            % Guess At Initial Conditions
[t,y] = ode45(@odes, tspan, ic);

figure
plot(t, y(:,1))
grid

% legend(compose('y(%d)',1:6), 'Location','best')


function dydt = odes(~,y)
x = y(1) ;
v = y(2) ;
temp = y(3) ;
temp_p = y(4) ;
%constants
c = .1 ; % t_c / t_h
a = .1 ; %alpha
p_i = 2 ; %initial p
temp_i = 1 ; %initial temp
temp_p_i = 1 ; %initial piston temp
friction_constant = .1 ;
% friction constant based on direction
f = -friction_constant*tanh(y(2)*10);
% % if y(2) < 0
% %     f = friction_constant ;
% % else
% %     f = -friction_constant ;
% % end
p = (temp ./ x) * p_i ;
dxdt = v ;
dvdt = p - 1 + f ;
dtempdt = -((c) * (temp - temp_p)) - ((1 / p_i ) * (2/3) * (p .* v)) ;
dtemp_pdt = (a * c * (temp - temp_p)) + ((a / p_i) * (2/3) * f * abs(v)) ;
dsgdt = ((3/2) * (1 ./ temp) .* dtempdt) + (v ./ x);
dspdt = (1/a) * (3/2) * ((1 ./ temp_p) .* dtemp_pdt) ;
dydt = [dxdt ; dvdt ; dtempdt ; dtemp_pdt ; dsgdt ; dspdt] ;
end

.

[Josh Ciesar]

I believe this solved the issue. Using the tanh function rather than an if statement is creative, I'll definitely use this again sometime. Thank you.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.