Trouble defining a set of parameters while using ode45 to solve a system of differential equations

I am trying to solve a system of differential equations where a piston is oscillating and colliding with a wall. In this system there are four differential equations:

p = (temp ./ x) * x_i ;

dxdt = v ; %x

dvdt = p ; %v

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

dtemp_pdt = a * c * (temp - temp_p) ; %piston temp

When I solve this system, it shows a piston with a single oscillation then it goes off into infinity, which is fine. Except I need to code the piston colliding with a wall. I drew a picture below to illustrate what I am talking about

When the piston collides with the right wal (x = 1), I need to code the following

v = - e * v ;

temp_p = temp_p + a * (2/3) * ((1 - (e^2))/2 * x_i) * (v.^2) ;

Where e, a ,x_i are constants. Only at x=1 these equations are true and they are supposed to define the new vleocity and piston temperature as the piston oscillates back to the left, and goes back to the first equations above. I am not sure how to include these equations since I am mixing non-differential equations with differential equations. How can I code this? The rest of my code is below.

%constants
c = .1 ; % t_c / t_h
a = .1 ; %alpha
p_i = 1 ; %initial p
temp_i = 1 ; %initial temp
temp_p_i = 1 ; %initial piston temp
x_i = .5 ;
e = .1 ;

y0 = [1 0 1 1] ;
t0 = 0 ; %initial t
tf = 100 ; %final t
tspan = [t0 tf] ;

[t,y] = ode45(@odes,tspan,y0) ;

x = y(:,1) ;
v = y(:,2) ;
temp = y(:,3) ;
temp_p = y(:,4) ;

figure(1)
plot(t,x)

figure(2)
plot(t,v)

figure(3)
plot(t,temp)

figure(4)
plot(t,temp_p)

%system of 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 = 1 ; %initial p
temp_i = 1 ; %initial temp
temp_p_i = 1 ; %initial piston temp
x_i = .5 ; %initial x
e = .1 ;

p = (temp ./ x) * x_i ;
dxdt = v ;
dvdt = p ;
dtempdt = -(c * (temp - temp_p)) - ((1 / x_i) * (2/3) * (p .* v)) ;
dtemp_pdt = a * c * (temp - temp_p) ;

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

end

Answer 1

[Torsten]

This page should help you:

https://de.mathworks.com/help/matlab/math/ode-event-location.html

[Josh Ciesar]

I have figure out how to stop integration when the event is located. But how do I restart integration with the new parameters?

[Torsten]

ODE45 returns control to the calling program after an event is located. Thereafter, you will have to call ODE45 again with modified initial conditions / equations / time span of integration. But I think this is all demonstrated in the ballode example, isn't it ?

[Josh Ciesar]

I managed to figure out how to do this. I used the link you attached as well as this answer here: https://www.mathworks.com/matlabcentral/answers/426186-confusion-on-event-function

Although, now I have another problem. In my case, there is a wall where the piston cannot go past. My code evaluates everything fine until the system approaches its final end point and then it looks like the wall disappears. I attached a graph below for reference.

This does not make sense since the velocity should become 0, and no more oscillations should occur. I realize that I could just set the xlim to where this error occurs, but I would like to continue the integration of other variables. Is this possible? How could I do this? All my code is below.

%constants

c = .1 ; % t_c / t_h

a = .1 ; %alpha

p_i = 1 ; %initial p

v_i = 0 ;

temp_i = 1 ; %initial temp

temp_p_i = 1 ; %initial piston temp

x_i = .5 ;

e = .9 ;

tf = 100 ;

y0 = [.5 v_i temp_i temp_p_i 0 0] ;

options = odeset('Events',@wallevent) ;

t = 0 ;

y = y0 ;

t_all = t ;

y_all = y ;

t_track = 0 ;

te = 0 ;

te_all = te ;

while t(end) < tf

[t,y,te,ye,ie] = ode45(@odes,[t_all(end) tf],y0,options) ;

t_all = [t_all ; t] ;

y_all = [y_all ; y] ;

te_all = [te_all ; te] ;

y0 = y(end,:) ;

y0(4) = y(end,4) + (a * (2/3) * ((1-(e^2))/(2 * x_i)) * (y(end,2)) .^2) ;

y0(2) = -1 * e * y(end,2) ;

end

x = y_all(:,1) ;

v = y_all(:,2) ;

temp = y_all(:,3) ;

temp_p = y_all(:,4) ;

s_g = y_all(:,5) ;

s_p = y_all(:,6) ;

p = (temp ./ x) * x_i ;

w_gas = - ((p .* v) / x_i) ;

w_p = (v / x_i) .* p ;

function dydt = odes(~,y)

x = y(1) ;

v = y(2) ;

temp = y(3) ;

temp_p = y(4) ;

s_g = y(5) ;

s_p = y(6) ;

%constants

c = .1 ; % t_c / t_h

a = .1 ; %alpha

x_i = 1 ; %initial x

e = .9 ;

p = (temp ./ x) * x_i ;

dxdt = v ;

dvdt = p ;

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

dtemp_pdt = a * c * (temp - temp_p) ;

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

function [position,isterminal,direction] = wallevent(t,y)

position = (y(1) - 1) ;

isterminal = 1 ;

direction = 1 ;

end

[Torsten]

The distance between the events when y(1) = 1 become smaller and smaller. Once ODE45 misses one such event, the reflection at the wall doesn't happen and the system develops as if no wall were present. I think at a certain time instant, you will have to change your equations to catch the phase after the velocity has settled to 0.

%constants
c = .1 ; % t_c / t_h
a = .1 ; %alpha
p_i = 1 ; %initial p
v_i = 0 ;
temp_i = 1 ; %initial temp
temp_p_i = 1 ; %initial piston temp
x_i = .5 ;
e = .9 ;
tf = 100 ;
y0 = [.5 v_i temp_i temp_p_i 0 0] ;
options = odeset('Events',@wallevent) ;
t = 0 ;
y = y0 ;
t_all = t ;
y_all = y ;
t_track = 0 ;
te = 0 ;
te_all = te ;
while t(end) < tf
    [t,y,te,ye,ie] = ode45(@odes,[t_all(end) tf],y0,options) ;
    t_all = [t_all ; t] ;
    y_all = [y_all ; y] ;
    te_all = [te_all ; te] ;
    y0 = y(end,:) ;
    y0(4) = y(end,4) + (a * (2/3) * ((1-(e^2))/(2 * x_i)) * (y(end,2)) .^2) ;
    y0(2) = -1 * e * y(end,2) ;
end
x = y_all(:,1) ;
v = y_all(:,2) ;
temp = y_all(:,3) ;
temp_p = y_all(:,4) ;
s_g = y_all(:,5) ;
s_p = y_all(:,6) ;
p = (temp ./ x) * x_i ;
w_gas = - ((p .* v) / x_i) ;
w_p = (v / x_i) .* p ;
figure(1)
plot(t_all(t_all<=20.5),y_all(t_all<=20.5,1))

figure(2)
plot(t_all(t_all<=20.5),y_all(t_all<=20.5,2))

function dydt = odes(~,y)
x = y(1) ;
v = y(2) ;
temp = y(3) ;
temp_p = y(4) ;
s_g = y(5) ;
s_p = y(6) ;
%constants
c = .1 ; % t_c / t_h
a = .1 ; %alpha
x_i = 1 ; %initial x
e = .9 ;
p = (temp ./ x) * x_i ;
dxdt = v ;
dvdt = p ;
dtempdt = -(c * (temp - temp_p)) - ((1 / x_i) * (2/3) * (p .* v)) ;
dtemp_pdt = a * c * (temp - temp_p) ;
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
function [position,isterminal,direction] = wallevent(t,y)
position = (y(1) - 1) ;
isterminal = 1 ;
direction = 1 ;
end