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