Is it possible to solve a nonlinear system with signum function using ODE45?
6 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
回答 (5 件)
Jan
2012 年 7 月 31 日
編集済み: 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.
Star Strider
2012 年 7 月 28 日
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 件のコメント
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!
Star Strider
2012 年 7 月 29 日
編集済み: Star Strider
2012 年 7 月 29 日
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 件のコメント
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.
shegaw mamo
2019 年 8 月 16 日
how can I get matlab command syntax or code for mathematical modeling of cancer treatment by radiotherapy
2 件のコメント
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.
参考
カテゴリ
Help Center および File Exchange で Symbolic Math Toolbox についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!