how to solve non linear differential equation with sign function

I have to solve the equation of the motion of a rigid block subjected to a horizontal base excitation ag. The equation of the motion is:

a + nk*g*sign(v) = - ag

a=acceleration of the block, nk=friction coefficient, g=gravity acceleration, v=velocity of the block, ag=horizontal base excitation

I am trying to solve it by ODE45 solver but because of the sign function it doesn't work: I have wrong results. I would like to know how to solve this kind of equation on Matlab. Thank you so much!

FUNCTION:

function xdot = sliding(a,x)

n=size(a);

nk = 0.3; % friction coefficient

g = 9.810; % Gravitational constant in m/sec^2

for i=1:n(1)

      xdot = [x(2);

              -a(i,1)-nk*g];
  end

SCRIPT:

a=0:0.1:4; %acceleration

t=0:0.1:4; %time

a=[a' t'];

x0 = [0; 0]; %initial condition

[t,x] = ode45('sliding',a(:,2),x0);

plot(t,a(:,1),'-');
xlabel('time'); ylabel('a(t)'); title('a (t)');
figure;
plot(t,x(:,1),'-'); xlabel('time'); ylabel('y_{1}(t)'); title('x (t)');
figure;
plot(t,x(:,2),'-'); xlabel('time'); ylabel('y_{2}(t)'); title('d x / dt (t)');

Answer 1

[Star Strider]

I knew I’d seen this problem before! See if this answer ( Is it possible to solve a nonlinear system with signum function using ODE45? ) applies to your problem.

[sara grossi]

Thank you so much for your reply! I saw that answer and I followed it for my script. Unfortunately it doesn't solve my problem. I obtain not realistic value and trend of the displacement and the velocity. In the last post it is said that integrating sign function with ODE45 gives errors and it is proposed an approach to solve the problem but I don't understand how implement it! Thank you so much!

[Star Strider]

There’s some missing information. I need to know how the differential equation changes when the sign changes. See Coulomb ODE with GUI for an example. (It has three conditions depending on what the sign() evaluates to.)

Without that information, it is only possible to integrate it to the first sign change, if I include an odeset options structure that stops it at the sign change. (This is correct.) The problem is that there is nowhere to go from there.

I got this to work, but it takes a while for it to integrate:

nk = 0.3;   % friction coefficient
g  = 9.810; % Gravitational constant in m/sec^2

x0 = [1;  0]*1E-5;

tstart = 0;
tfinal = 4;

sliding = @(t,x) [x(2);  -x(1)-nk*g*sign(x(2))];

[t,x] = ode45(sliding,[tstart tfinal],x0);

I’m not sure it’s correct though, because as written, it’s integrating across the discontinuities. It should restart after each sign() change with new initial conditions and a new value for the Coulomb friction force.

[sara grossi]

Thank you so much for your reply!

I am sorry: I haven't explained the whole problem. The body rests until the acceleration of the ground ag is bigger than a costant value c: ag>c then the equation of the motion is:

a(t) = - ag - nk*g*sign(v(t))

sign(v)=1 if v>0, sign(v)=-1 if v<0. If v=0 the body rests and the problem starts from the biginnig: ag>c the body will starts move otherwise it will rest. In your code you put x(1) in the sliding function instead ag that is an external input. (I am sorry in my script it was not clear that aspect!). Following the example "Coulomb ODE with GUI" that you suggested I am trying to change my code. (for now I have given a linear ground acceleration.)

SCRIPT:

ag=0:0.1:10; %acceleration of the ground
t=0:0.1:10; %time in which the problem is solved
n=size(ag);
y0 = [0; 0]; %initial condition
y = zeros(2,n(2)+1); %Initialize the solution vector
y(:,1) = y0; %Set the initial conditions   
for k = 1:n(2)
    if abs(ag(:,k))<=0.39;
            y(:,k+1)=y0;
    else [T y(:,k+1)] = ode45(@(t,x) cnsegno(t,y(:,k+1),ag(:,k)),t,y(:,k));
    end
end

FUNCTION:

function dxdt = cnsegno(t,xs,a)
nk = 3; % coeff att
g = 9.810; % Gravitational constant in mm/sec^2

dxdt = [xs(2);
        -ag-sign(xs(2))*nk*g];

I have give the function as an array because the input motion, ag, is a coefficient that depend on the time. Right now it still doesn't work and I don't know how to fix it. Could you help me? In addition I don't know how to stop the analysis when v=0 and then restart the analysis when the ground acceleration became bigger than c again. Thank you so much for your help. I really appreciate it!!

[Star Strider]

Thank you for the additional information. I will review your posts and see if I can come up with a coherent solution.

So the system is at rest until a*g > c. Then if v(t) = 0, that restarts the integration. (That may be the information I need for my ‘event’ function.) I could not figure out what the function was supposed to do at the discontinuities. I’ll work on an odeset ‘event’ options structure and function, and get back to you.

I do not understand how v(t) can change signs and not trigger a zero-crossing event, though. I do not see a value for c, so I will wait for it before I continue this.

I do not understand your a vector though. I thought acceleration, x"(t) = a was supposed to be a value the ODE calculated and returned.

I apologise, but I am a bit lost.

[sara grossi]

I'm so sorry if I haven't explained properly!

My a vector is ag: the acceleration of the ground. (I have just corrected my previous comment). I have to find the acceleration of the body that is expressed by the equantion: a(t) = - ag - nk*g*sign(v(t)).

c=3,9

v(t) can't change signs without trigger a zero-crossing event. And when v(t) become zero the problem starts from the biginnig: if ag<c the body rests, otherwise the body starts move again and the ode has to be solved again.

Do you know why the ode takes so long to work?

Thank you so much for yuor help!

[Star Strider]

I thought we needed three functions of time:

• position, x(t)
• velocity, x’(t)
• acceleration x”(t)

That makes a(t) = x”(t), correct? That is how I am used to thinking of problems like this. I think part of my problem is understanding the background of this problem.

I still do not understand ‘ ag: the acceleration of the ground. ’ I also do not understand why it is a vector.

I apologise, but I am still lost.

[sara grossi]

I apologize again for not being clear. I rephrase the problem:

The motion is initiated from rest once ag > 3.9. (if ag < 3.9 the body rests.)

Once sliding, the governing equation of motion may be expressed as x''(t) = -ag(t) - S(x'(t))*nk*g in which S(x'(t)) is defined through the signum function: S(x'(t))=1 if x'(t)>0, S(x'(t))=-1 if x'(t)<0 The motion continues until the relative velocity of the body equals zero: x'=0 At that time, the body momentarily comes to rest relative to the ground and ag > or < 3.9 is evaluated to determine if it remains at rest, or starts to slide again.

ag is the input that makes move the body. It is the external acceleration given by the ground at the base of the body. ag is a vector because it is time dependent ag(t). When my code will be ready I will give it different ag. Rigth now I am trying to solve for a linear ground acceleration.

x''(t) = acceleration of the body relative to the ground; x'(t) = velocity of the body relative to the ground; nk= friction coefficient; g=gravity acceleration.

So there are two acceleration: the acceleration of the ground ag(t) and the acceleration of the body x''(t). I need to solve the equation to have x'(t), x(t).

Thank you for your patience!