ode45 with discrete data

My problem consists of a differential equation in which, the damping coefficient of my system is not linear, it depends on the displacement

where, on the x axis is the displacement of the body and on the y axis is the damping coefficient for these positions.

I did the code in matlab and he did the calculation, however, I'm not sure if he did it correctly, since, in the damping coefficient curve, there is no coefficient for all positions.

An important point is, that I used the slipne function to get the equation, so, I wanted your help to understand a little more about the equation

And my question is more for, I know that the displacement of the system will be within the range that is in the graph of the damping coefficient, but does matlab understand that? why, although I haven't provided all the damping coefficient values, will the matlab integrate only in the range that I provided?

The other point is, I plotted the spline function and, from very high displacements, for my system, the damping coefficient starts to have negative values, so, I want to understand how the matlab can read this function and if it solves my problem in the range provided.

function testando

k1 = 116.0665; % N/m
m1 = 0.06; % kg
w = 25.1327; % rad/s
F1 = 0.01;% N
f1 = @(time) F1.*sin(w*time);
     
x = xlsread('emf.xlsx','B2:B62');
FT = xlsread('emf.xlsx','J2:J62');

x0 = 0; v0 = 0;

IC = [x0,v0];

% Time

t0 = 0; tf = 10;

tspan = [t0,tf];

%sdot = g(t,s)
sdot1 = @(t,s) [s(2);
    (f1(t) - spline(x,FT,s(1)).*s(2) - k1.*s(1))./m1];
% Integração numerica 

[time,state_values] = ode45(sdot1,tspan,IC);
displacement = state_values(:,1);

% plots

plot(time,displacement),xlabel('time(s)'),ylabel('displacement(m)')
title('Displacement(m) vs. t')

end

Answer 1

[Ameer Hamza]

For your first question, yes, the spline function only uses the displacement value from vector 'x' the coefficient value from vector 'FT'. As long as the displacement value s(1) lies in the range of 'x', the MATLAB will use values in the range you provided.

For the second question, you are correct. The spline function can overshoot from the actual values. I recommend using interp1, which uses linear interpolation by default. The interpolated value will exactly look like the figure you attached.

sdot1 = @(t,s) [s(2);
    (f1(t) - interp1(x,FT,s(1)).*s(2) - k1.*s(1))./m1];

[José Anchieta de Jesus Filho]

Thanks for your response

So, let me see if I understand, I gave the matlab an interval in which, in that interval, the damping coefficient will assume different values for each displacement. I know that the solution, speaking of the real model of the system, that it will be within that range, so MATLAB will use only the range that I provided, so I don't have to worry about the inaccuracies of cubic interpolation for large displacements, correct ?

Another point is, I tried to use interp1, but it did not return a solution, there was no graph. He managed to solve it using, for example, 'spline', 'pchip' and 'cubic'

sdot1 = @(t,s) [s(2);
    (f1(t) - interp1(x,FT,s(1),'cubic').*s(2) - k1.*s(1))./m1];

[José Anchieta de Jesus Filho]

If, for example, I need to place an "if and else" condition to vary the value of that coefficient.

Ex: if s (1) is between -0.0302 and 0.0298, the value of the coefficient is provided in the function, if it is not, the value must be 0. How could I perform this logical operation in MATLAB

[Ameer Hamza]

1. Yes, If the value you provided is within the range of the input displacement, then the cubic interpolation will not have much effect. You can see the effect on the interpolation by making a graph using this code

x % your matrix x
FT % your matrix Ft
xq = linspace(min(x), max(x), 1000);
FTq = spline(x,FT,xq);
plot(x,FT,'+',xq,FTq);

This will show how the spline function interpolates in-between values.

2. This is strange behavior with interp1 with linear interpolation. Can you use the following code to see if the interpolation is correct

x % your matrix x
FT % your matrix Ft
xq = linspace(min(x), max(x), 1000);
FTq = interp1(x,FT,xq);
plot(x,FT,'+',xq,FTq);

3. For that you will need to write a seperate function

function dsdt = odeFun(t,s,x,Ft)
    if -0.0302 < s(1) && s(1) < 0.0298
        c = spline(x,FT,s(1));
    else
        c = 0;
    end
    
    sdot1 = @(t,s) [s(2);
        (f1(t) - spline(x,FT,s(1)).*s(2) - k1.*s(1))./m1];
end

and then call ode45 like this

[time,state_values] = ode45(@(t,s) sdot1(t,s,x,Ft),tspan,IC);

[José Anchieta de Jesus Filho]

Thank you very much, your codes helped to understand a little more about ode45.

I tested the first function I did, without any conditions, and then I tested the function using the "if and else" conditions. The result was the same.

Really, matlab understood that only the interval I provided was necessary.

Thank you so much for this.

function testando

k1 = 116.0665; % N/m
m1 = 0.06; % kg
w = 25.1327; % rad/s
F1 = 0.01;% N
f1 = @(time) F1.*sin(w*time);
     
x = xlsread('emf.xlsx','B2:B62');
FT = xlsread('emf.xlsx','K2:K62');

x0 = 0; v0 = 0;

IC = [x0,v0];

% Time

t0 = 0; tf = 10;

tspan = [t0,tf];

%sdot = g(t,s)
sdot1 = @(t,s) [s(2);
    (f1(t) - interp1(x,FT,s(1),'spline').*s(2) - k1.*s(1))./m1];
% Integração numerica 

[time,state_values] = ode45(sdot1,tspan,IC);
displacement = state_values(:,1);

% plots
figure(2)
plot(time,displacement),xlabel('time(s)'),ylabel('displacement(m)')
title('Displacement(m) vs. t')

end

The code using the conditions

function sdot = odeFun(t,s,x,Ft)
k1 = 116.0665; % N/m
m1 = 0.06; % kg
w = 25.1327; % rad/s
F1 = 0.01;% N
f1 = @(time) F1.*sin(w*time);
    if -0.0302 < s(1) && s(1) < 0.0298
        c = spline(x,Ft,s(1));
    else
        c = 0;
    end
    
    sdot(1) = s(2);
    sdot(2) = (f1(t) - c.*s(2) - k1.*s(1))./m1;
    sdot = sdot';
        
end

[Ameer Hamza]

I am glad to be of help.