-2

How can the non-smooth functions be integrated using Matlab? Can the Cash-Karp Method be used? The Dormand-Prince method seems to give an error while integrating a nonsmooth solution.

I am using MATLAB ODE45 for the integration of my ordinary differential equation.

enter image description here

$\{q\}$ is a Column Vector of length 14. The excitation force on the right-hand side varies in the following form

enter image description here

where $\{h_p\}$, $\{g_p\}$, $\{h_g\}$, and $\{g_g\}$ are column vectors of length 6

I call it nonsmooth because $\{F_\text{excit}\}$ varies with $F_\text{mf}$.

The system of equations can be solved when $F_\text{mf}=0$. When $F_\text{mf}$ is nonzero, I am getting a warning which actually indicates that integration cannot be performed.

Warning/Error Message in Matlab

Warning: Failure at t=1.821477e+00.  Unable to meet
integration tolerances without reducing the step 
size below the smallest value
allowed (3.552714e-15) at time t.

I did try ODE15S but the error appeared much earlier than ODE45. Both gave me the same error

I use square brackets to represent matrices and curly braces to represent vectors.

Anton Menshov
  • 8,672
  • 7
  • 38
  • 94
Srikumar Gopalakrishnan
  • 3
  • 3
  • 1
    Can you be more explicit about what exactly is not smooth? – Wrzlprmft Nov 13 '16 at 07:17
  • 6
    Please post the equations you're trying to solve and (exactly) the error you get. – David Ketcheson Nov 13 '16 at 13:58
  • https://scicomp.stackexchange.com/questions/24998/which-runge-kutta-method-is-more-accurate-dormand-prince-or-cash-karp – Mathews24 Nov 13 '16 at 22:56
  • 3
    What do the braces and brackets in your notation mean? If they are just normal brackets, why does $F_\text{excit}$ just have four components? — I call it nonsmooth because Fexcit varies with Fmf – I fail to see how this makes anything nonsmooth. – Wrzlprmft Nov 14 '16 at 17:58

2 Answers2

3

You really have not provided enough information to diagnose the cause of your problem (e.g. I don't know if the dimension of your q-vector is two or two million) but guessing, based on your error message, I doubt it has anything to do with "smoothness" of the forcing function.

Generally, the way to debug these MATLAB ODE solvers when you get the error message you show is to repeat the analysis, stopping just before the failure time (e.g. 1.8 seconds). You want to save the solution at enough time steps so you can plot some of the key dependent variables as a function of time. I think it is very likely that one or more of these is going to plus or minus infinity near this final time.

The reason for this is most likely an error in your formulation or input to ode45.

A second possibility is that your system is inherently unstable.

A third possibility is that your system of ODE is stiff and that is causing problems for ode45. ode45 is not designed to solve stiff ODE systems. Changing to an ode solver with a different Runge-Kutta pair or some other explicit solver will not solve this problem. The MATLAB solver-of-choice for stiff systems of the form you show is ode15s.

Bill Greene
  • 6,064
  • 1
  • 16
  • 24
  • Can you explain what is smoothness of the forcing function ?, I was referring to a paper on ACM by Cash & Karp. If you have any inputs on the references to read that will be of great help – Srikumar Gopalakrishnan Nov 14 '16 at 20:05
  • You were the one who mentioned smoothness of the forcing function. I simply said I doubted if that was the cause of the premature termination of ode45. I'm not that surprised you saw the same behavior with ode15s; stiffness was my third suggestion as to the cause of your problem. The fact that you see the same behavior with ode15s is good evidence that you won't see different behavior simply by changing the Runge-Kutta form slightly. – Bill Greene Nov 14 '16 at 20:27
  • But the problem looks is due to the second possibility you stated. Is there any reference which explains about stability criteria for an iterative system. To ensure the stability the spectral radius of the system (max eigen value of the iterative sytem) should be less than one. How to determine the spectral radius of the forced vibration system which has damping in it ?, Any reference for this will be highly appreciated. – Srikumar Gopalakrishnan Nov 15 '16 at 04:07
  • I don't know what you mean by "iterative system." Maybe you are confusing stability of a mechanical system with that of a numerical procedure. Most books on dynamical systems discuss stability. – Bill Greene Nov 15 '16 at 12:41
  • If all of the matrices on the LHS of your equation are positive definite, your system will be stable. It is likely that it is the particular form of the excitation on the RHS that is causing the response to become very large. – Bill Greene Nov 15 '16 at 15:03
0

I am afraid to have understood the nonsmoothness of your problem correctly, but since it looks like a typical equation of motion in mechanics, so it may be similar to friction and impacts. For these kind of ODEs exist two classes of numerical solution: event-driven and time-stepping, whereas standard methods assuming smoothness fail. For the event-driven there is an example (bouncing ball) in the MatLab documentation, for time-stepping (aka event-collecting), I give a simple example (I am not aware of a built-in solver in MatLab)

% Implementation of a time−stepping scheme [Moreau1988]
% for a bouncing ball (1D), ground modeled by unilateral constraint
% y>=0 and coefficient of restitution dy/dt plus=−rc ∗ dydt minus
clear; clc; close all;
N=1000; h=0.01 ; T=N∗h ; t=0:h:T; % time discretization
rc=0.9; % coefficient of restitution
m=1; g=10;
 % mass and gravitational acceleration
y=zeros(N+1, 1); v=zeros (N+1, 1) ; % position and velocity arrays
R=zeros(N, 1); % reaction impulse, value at t=0 depends on past
F=zeros(N, 1); % external impulse, value at t=0 does not enter
% i.c . and forcing for bouncing −−> rest "contact closing" 
q0=0.5; v0=0; % initial conditions
f=@(t) −m∗g; % external force (continuous−time)
% a sum of a geometric series gets a physical meaning [Zeno's paradoxes]
zenotime=sqrt(8∗q0/g)/(1−rc)−sqrt(2∗q0/g); %assuming v0=0 
% i.c. and forcing for resting −−> flight "contact opening" 
%q0=0.0; v0=0; % initial conditions
%f=@(t) (−1.2∗sin(t∗2∗pi/T−1)∗m∗g; % external force (continuous−time)
y(1)=q0; v(1)=v0;
for n=1:N
    % Euler−forward for y(n+1)=y(n)+v(n)∗h; position
    % Euler−backward for velocity
    F(n)=h∗f(t(n+1));% integrated force at t=t(n+1)
    v(n+1)=v(n)+F(n)/m; % free flight, otherwise overwritten by LCP
    % LCP for contact, manually resolved
    if (y(n+1)<=0)
        % prevent penetration
        dv=−(1+rc)∗v(n); % (potential) velocity jump
        if (m∗dv)>F(n) % ground can only push, but not pull
            v(n+1)=v(n)+dv; % impact law
            R(n)=m∗dv−F(n); % reaction force from equation of motion
        end % else nothing to do
    end % else nothing to do
end
figure;
yyaxis left; plot(t, y, [zenotime, zenotime], [−0.1, 0.6], 'k−−') ;
xlabel('time [s]'); ylabel('position [m]'); ylim([−0.1 0.6])
yyaxis right; plot(t(2:end), R); % R at t=0 is not determined
ylabel ('reaction impulse [Ns]'); ylim ([−1 6] );

For friction I may provide a similar example, too.

Dominik Kern
  • 86
  • 6