Is it possible to solve a nonlinear system with signum function using ODE45?

Hi!

I'm working with friction, my system is a SDOF with a mass, stiffness, friction and a harmonic function applied, then the equation of motion can be expressed by: ma+kx=f*sin(wt) but depending on the sign of the relative velocity, the Coulomb friction term should be added: Fd=ff*sign(velocity).

I would like to know if this differential equation is possible to be calculated using the solver ode45 or if it is necessary to check the relative velocity at each time step to use the correct equation.

Thank you

Laura

ode45 solution:

xp=zeros(2,1);
xp(1,1)=x(2,1);
xp(2,1)= (- K * x(1,1) - Fd*sign(x(2,1)) + (Fi * sin(wi*t)))/M;

Answer 1

[Jan]

No, do not use ODE45 to integrate functions, which are not continuously differentiable. Otherwise you break the specifications and the stepsize control will jump to its limits until the local discretization error is dominated by rounding errors. Of course you can get a value out of this integration, but it can be dominated by random rounding effects.

A scientifically clear approach is to stop the integration at every discontinuity in the integrated function and use the current positions (plus modifications by the jump on demand) to restart the integration. The event functions and an additional outer loop around ODE45 allow this by a few additional lines of code.

Although it is a mistake from the point of proper numerical calculations, you find "results" obtained by the integration of discontinuous functions in a lot of papers, diploma theses and dissertation (apart from the field of numerics). Unfortunately this is not mentioned in the help or doc of ODE45, but the documentation of Matlab cannot (and should not) replace a participation in classes for numerics.

[laura.gcas]

Thank you. I didn't know that so well. I'm quite new using matlab and numerical analysis.

Now I have created a similar code for this simple case using Newmark integration and using two conditions: when the friction is considered and when it is not, therefore I haven't used the solver ODE45.

However, in this case the results using ODE45 and Newmark are almost the same. As you said they might have some errors but I guess these errors may become worse if I want to increase the dofs of the system.

Thanks again for your comments

Answer 2

[Star Strider]

Put it in a function and give it a go! The statement looks as if it's correctly coded.

Let us know if you have problems. If you do, post your code and a short description of what doesn't work, any error messages you get, what lines they refer to (list the lines specifically because they're not numbered here), and an example of the output you get. Remember to format it as:

Code

Also include your initial conditions, time span, and values for K, Fc, Fi, and M so we can simulate it if necessary.

[Star Strider]

I coded it as an anonymous function and it works fine! I'd still like to know the values you're using for the constants. (I made some up to do the simulation.)

Have fun with it!

Answer 3

[laura.gcas]

Hi!! thanks for your reply

I have some problems using my values (I took them from an example), the code is the following:

M=1; %kg
k=1000; %N/m
Fd=1; %N
Fi=15; %N
wi=13; %rad/s

%%Initial values
x10=0;
x20=0;
x0=[x10 x20]';

% Simulation time
dt=0.01;
n=1000; 
tspan=linspace(0,n*dt,n);

% ODE45
pars=[M k m Fd kd c1 c2 Fi wi];
options=odeset('RelTol',1e-6);
[t,x]=ode15s(@friction_lastfunc,tspan,x0,options,pars);

function xp=friction_lastfunc(t,x,pars)
%%%model parameters
M=pars(1);
k=pars(2);
Fd=pars(3);
Fi=pars(4);
wi=pars(5);

KK = k;
MM = M;

xp=zeros(2,1); 
xp(1,1)=x(2,1);
xp(2,1)= (- KK * x(1,1) - Fd*sign(x(2,1)) + (Fi * sin(wi*t)))/M;

Does it work for you?

Thank you in advance

Laura

Answer 4

[Star Strider]

Your code and values work perfectly for me in my implementation of them (I didn't change your code), and produces a really interesting plot:

        M=1; %kg 
        K=1000; %N/m 
        Fd=1; %N 
        Fi=15; %N 
        wi=13; %rad/s

        CulombFrctn = @(t,x) ([ x(2,1); (- K * x(1,1) - Fd*sign(x(2,1)) + (Fi * sin(wi*t)))/M]);

        t0 = clock
        [t x] = ode45(CulombFrctn, [0 : 0.003 : 2], [0; 0]);
        t1 = clock;
        runtime = etime(t1,t0);

        xstats = [mean(x);  median(x);  std(x);  min(x);  max(x)]
        limsy = [xstats(2,2)-3*xstats(3,2)  xstats(2,2)+3*xstats(3,2)]

        X1Scl = fix(log10(xstats(3,2)/xstats(3,1)));

        figure(128)
        plot(t, x(:,1)*10^X1Scl, '-b', 'LineWidth',2)
        hold on
        plot(t, x(:,2), '-r', 'LineWidth',1)
        hold off
        legend(sprintf('x_{1} x 10^{%d}',X1Scl), 'x_{2}', 1)
        axis([xlim  xstats(2,2)-3*xstats(3,2)  xstats(2,2)+3*xstats(3,2)])
        grid

It took less than a second for ‘ode45’ to solve.

[laura.gcas]

Hi!! Thank you very much for your reply! It works. Before I created a function file where I had the ode and it was called from the main file, is that possible to cause problems? Because it is strange that adding the same function directly to the main file everything works fine...

Can I ask you why you use some statistical parameters in the code?: xstats = [mean(x); median(x); std(x); min(x); max(x)] limsy = [xstats(2,2)-3*xstats(3,2) xstats(2,2)+3*xstats(3,2)]

The plots look really nice to me :) Thanks a lot!

Laura

[Star Strider]

I calculated the statistics to set the limits of the plot because early attempts with the parameters I guessed at created some extreme values, and so that x(:,1) and x(:,2) were both large enough to be easily visible on the plot.

Answer 5

[shegaw mamo]

how can I get matlab command syntax or code for mathematical modeling of cancer treatment by radiotherapy

[Walter Roberson]

https://www.mathworks.com/company/mathworks-stories/cancer-therapy-targets-cancer-in-real-time.html

https://www.mathworks.com/matlabcentral/answers/472438-matlab-saber-spatial-and-biological-evaluation-for-radiotherapy

https://www.mathworks.com/academia/books/clinical-radiotherapy-physics-with-matlab-dvorak.html

[Jan]

@hegaw mamo: Please open a new thread to ask a new question. Attaching a question in the section for answers of another thread is less useful. Now you cannot accept answers as solution and otehr readers will not find your question.