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

Question

laura.gcas 2012-7-28

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/44752-is-it-possible-to-solve-a-nonlinear-system-with-signum-function-using-ode45

评论： Jan 2019-8-19

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;

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

请先登录，再进行评论。

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

Answer 1

Jan 2012-7-31

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/44752-is-it-possible-to-solve-a-nonlinear-system-with-signum-function-using-ode45#answer_55041

编辑：Jan 2012-7-31

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

laura.gcas 2012-7-31

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 2012-7-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/44752-is-it-possible-to-solve-a-nonlinear-system-with-signum-function-using-ode45#answer_54830

在 MATLAB Online 中打开

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Star Strider 2012-7-29

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 2012-7-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/44752-is-it-possible-to-solve-a-nonlinear-system-with-signum-function-using-ode45#answer_54850

编辑：Walter Roberson 2012-7-31

在 MATLAB Online 中打开

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

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

请先登录，再进行评论。

Answer 4

Star Strider 2012-7-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/44752-is-it-possible-to-solve-a-nonlinear-system-with-signum-function-using-ode45#answer_54863

编辑：Star Strider 2012-7-29

在 MATLAB Online 中打开

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.

2 个评论
显示无隐藏无

laura.gcas 2012-7-31

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 2012-7-31

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 2019-8-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/44752-is-it-possible-to-solve-a-nonlinear-system-with-signum-function-using-ode45#answer_387806

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

2 个评论
显示无隐藏无

Walter Roberson 2019-8-16

编辑：Walter Roberson 2019-8-16

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 2019-8-19

@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.

请先登录，再进行评论。