ODE45: Solving a system of ODEs, except one of the equations changes based on certain parameters. How can I code this?

I am trying to solve a problem regarding a piston with friction.

There are 5 equations regarding this system:

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 * v) ;

The 5th equation defines "p" or the pressure.

When the piston is moving in the rightward direction, the value of ' f ' in the function needs to be negative. When the piston is moving leftward, the value ' f ' needs to be positive. I have solved similar problems regarding this piston without friction, so I am aware of how to solve the system using ode45. However, I am not sure how I can code the changing value of ' f ', since the equation changes depending on the physical parameters. How could I code this?

Answer 1

[William Rose]

@Josh Ciesar, Is pressure, p, a state variable, given by a first order differential equation, or is it an instantaneous function of the state variables (i.e. not defined by a differential equation)? Please specify the equation for p.

Is it correct that "piston moving leftward" means v is negative, and "piston moving rightward" means v is positive? Are you sure that is correct? It makes me think that this "friction" is a force that adds to velocity when velocity is one way, and it opposes movement when movement is the other way. This does not sound like normal friction, which opposes velocity in either direction.

Getting back to your question, you can write the function that returns the derivative so that f is positive when v is positive, and negative when v is negative (or vice versa).

Here is an example. In this example, I assume that p is not a part of the state vector and that p=x*v*temp. I know that the second part of that statement is not correct. I realize the first part may also b incorrect. If you provide the equation for p, we can adjust things.

function dydt = derivative(t,y)
%y(1:4)=[x;v;temp;temp_p]
a=1;
c=2;
p_i=3;
if y(2)<0
    f=-10;  %f<0 if v<0
else 
    f=+10;  %f>0 if v>0
end
dydt=zeros(4,1);
p=y(1)*y(2)*y(3);
dydt(1)=y(2);
dydt(2)=p-1+f;
dydt(3)=-c*(y(3)-y(4))-(1/p_i)*2*p*y(2)/3;
dydt(4)=a*c*(y(3)-y(4))+(a/p_i)*2*f*v/3;
end

Try it. Good luck.

Answer 2

[Sam Chak]

@Josh Ciesar

Assuming that the right direction of the piston is defined by , then this graph shows exactly what you described.

x = linspace(-1, 1, 1001);
m = 1;          % defines the negativity slope of f
f = - m*x;      % only you how negative the f is.
plot(x, f, 'linewidth', 1.5), grid on, xlabel('Piston direction'), ylabel('f')

Else, it can look like this:

k = 1;              % magnitude of the f
f = - k*sign(x); 
plot(x, f, 'linewidth', 1.5), grid on, xlabel('Piston direction'), ylabel('f')
ylim([-2 2])