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

Question

Josh Ciesar 2022-8-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1772100-trouble-defining-a-set-of-parameters-while-using-ode45-to-solve-a-system-of-differential-equations

编辑： Torsten 2022-8-5

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-8-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1772100-trouble-defining-a-set-of-parameters-while-using-ode45-to-solve-a-system-of-differential-equations#answer_1019320

This page should help you:

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Josh Ciesar 2022-8-3

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 2022-8-5

编辑：Torsten 2022-8-5

在 MATLAB Online 中打开

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

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论