1#
2#     This file is part of CasADi.
3#
4#     CasADi -- A symbolic framework for dynamic optimization.
5#     Copyright (C) 2010-2014 Joel Andersson, Joris Gillis, Moritz Diehl,
6#                             K.U. Leuven. All rights reserved.
7#     Copyright (C) 2011-2014 Greg Horn
8#
9#     CasADi is free software; you can redistribute it and/or
10#     modify it under the terms of the GNU Lesser General Public
11#     License as published by the Free Software Foundation; either
12#     version 3 of the License, or (at your option) any later version.
13#
14#     CasADi is distributed in the hope that it will be useful,
15#     but WITHOUT ANY WARRANTY; without even the implied warranty of
16#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17#     Lesser General Public License for more details.
18#
19#     You should have received a copy of the GNU Lesser General Public
20#     License along with CasADi; if not, write to the Free Software
21#     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22#
23#
24#! Integrator tolerances
25#! =====================
26from casadi import *
27from numpy import *
28from pylab import *
29
30x=SX.sym('x')
31dx=SX.sym('dx')
32states = vertcat(x,dx)
33
34dae={'x':states, 'ode':vertcat(dx,-x)}
35
36tend = 2*pi*3
37ts = linspace(0,tend,1000)
38
39tolerances = [-10,-5,-4,-3,-2,-1]
40
41figure()
42
43for tol in tolerances:
44  opts = {'reltol':10.0**tol, 'abstol':10.0**tol, 'grid':ts, 'output_t0':True}
45  F = integrator('F', 'cvodes', dae, opts)
46  res = F(x0=[1,0])
47
48  plot(ts,array(res['xf'])[0,:].T,label='tol = 1e%d' % tol)
49
50legend( loc='upper left')
51xlabel('Time [s]')
52ylabel('State x [-]')
53show()
54
55
56tolerances = logspace(-15,1,500)
57endresult=[]
58
59for tol in tolerances:
60  opts = {}
61  opts['reltol'] = tol
62  opts['abstol'] = tol
63  opts['tf'] = tend
64  F = integrator('F', 'cvodes', dae, opts)
65  res = F(x0=[1,0])
66  endresult.append(res['xf'][0])
67
68figure()
69loglog(tolerances,(array(endresult)-1),'b',label='Positive error')
70loglog(tolerances,-(array(endresult)-1),'r',label='Negative error')
71xlabel('Integrator relative tolerance')
72ylabel('Error at the end of integration time')
73legend(loc='upper left')
74show()
75