How to numerically test stability of an oscillation curve? [closed]

Question

Closed. This question is off-topic. It is not currently accepting answers.

---

Want to improve this question? Update the question so it's on-topic for Stack Overflow.

Closed 8 years ago.

Improve this question

I am now having an oscillation curve which is part of the solutions of a set of nonlinear ordinary differential equations. I am required to test the stability/convergence of this curve as time goes to infinite. How to do it with Matlab?

The figure looks like this:

Answer 1

I turned out using the following script, it works fine for me, but I am still wondering, is there any better way of predicting convergence at long time.

function err = stability_test(t, y)
% Given data of an oscillating curve y(t), tell whether the oscillation 
%   amplitude decrease or not by 
%   1. locating peaking points
%   2. linear fit peak points and see if the gradient is negative or not
%
% t, y must be of the same shape
% err = 0, non-ocillating
%     < 0, stable
%     > 0, unstable
nt = linspace(min(t), max(t), 500);
ny = interp1(t,y,nt,'spline');
ndy = gradient(ny,nt);
ndy2 = del2(ny,nt);
if(isempty(find(ndy<0, 1)) || isempty(find(ndy2>0, 1)))
    err = 0;
else
    ndt = nt(2) - nt(1);
    ii = find(abs(ndy)<abs(ndt*ndy2*2) & ndy2<0);

    if(isempty(ii))
        err = 0;
    else
        if(length(ii)==1)
            ii = [ii,length(ndy)];
        end
        ym = ny(ii);
        tm = nt(ii);
        p = polyfit(tm, ym,1);
        err = p(1);
    end
end

Answer 2

It has been eight years since I did anything like this, so take my answer with a grain of salt.

1. Solve the equations using step size S and also with step size S/2; if the results match (i.e. are within machine epsilon, or 10x machine epsilon, or however you're defining the word "match"), then you're good to go on truncation error
2. Solve the equations using standard floating point arithmetic, and also solve them with extended precision arithmetic (IIRC Matlab calls this Variable Precision Arithmetic; IEEE double precision arithmetic uses a 52-bit significand, so an 80-bit significand ought to be more than enough to reveal instability due to roundoff error); if the results match, then you're good to go on roundoff error