diff --git a/.codedocs b/.codedocs new file mode 100644 index 0000000..60c6fb2 --- /dev/null +++ b/.codedocs @@ -0,0 +1,193 @@ +# CodeDocs.xyz Configuration File +# +# Rename this example to '.codedocs' and put it in the root directory of your +# repository. This file is optional, documentation will still be generated +# without it using sensible defaults. + +#--------------------------------------------------------------------------- +# CodeDocs Configuration +#--------------------------------------------------------------------------- + +# Include the Doxygen configuration from another file. +# The file must be a relative path with respect to the root of the repository. +# If any of the options in this doxyfile include a path (ie, INPUT), these +# paths will be considered relative to the root of the repository, not the +# location of the DOXYFILE. + +DOXYFILE = Doxyfile + +# Specify external repository to link documentation with. +# This is similar to Doxygen's TAGFILES option, but will automatically link to +# tags of other repositories already using CodeDocs. List each repository to +# link with by giving its location in the form of owner/repository. +# For example: +# TAGLINKS = doxygen/doxygen CodeDocs/osg +# Note: these repositories must already be built on CodeDocs. + +TAGLINKS = + +#--------------------------------------------------------------------------- +# Doxygen Configuration +#--------------------------------------------------------------------------- + +# Doxygen configuration may also be placed in this file. +# Currently, the following Doxygen configuration options are available. Refer +# to http://doxygen.org/manual/config.html for detailed explanation of the +# options. To request support for more options, contact support@codedocs.xyz. +# +# ABBREVIATE_BRIEF = +# ALIASES = +# ALLEXTERNALS = +# ALLOW_UNICODE_NAMES = +# ALPHABETICAL_INDEX = +# ALWAYS_DETAILED_SEC = +# AUTOLINK_SUPPORT = +# BRIEF_MEMBER_DESC = +# BUILTIN_STL_SUPPORT = +# CALLER_GRAPH = +# CALL_GRAPH = +# CASE_SENSE_NAMES = +# CITE_BIB_FILES = +# CLASS_DIAGRAMS = +# CLASS_GRAPH = +# COLLABORATION_GRAPH = +# COLS_IN_ALPHA_INDEX = +# CPP_CLI_SUPPORT = +# DIAFILE_DIRS = +# DIRECTORY_GRAPH = +# DISABLE_INDEX = +# DISTRIBUTE_GROUP_DOC = +# DOTFILE_DIRS = +# DOT_FONTNAME = +# DOT_FONTSIZE = +# DOT_GRAPH_MAX_NODES = +# DOT_IMAGE_FORMAT = +# DOT_TRANSPARENT = +# DOXYFILE_ENCODING = +# ENABLED_SECTIONS = +# ENABLE_PREPROCESSING = +# ENUM_VALUES_PER_LINE = +# EXAMPLE_PATH = +# EXAMPLE_PATTERNS = +# EXAMPLE_RECURSIVE = +# EXCLUDE = +# EXCLUDE_PATTERNS = +# EXCLUDE_SYMBOLS = +# EXPAND_AS_DEFINED = +# EXPAND_ONLY_PREDEF = +# EXTENSION_MAPPING = +# EXTERNAL_GROUPS = +# EXTERNAL_PAGES = +# EXTRACT_ALL = +# EXTRACT_ANON_NSPACES = +# EXTRACT_LOCAL_CLASSES = +# EXTRACT_LOCAL_METHODS = +# EXTRACT_PACKAGE = +# EXTRACT_PRIVATE = +# EXTRACT_STATIC = +# EXT_LINKS_IN_WINDOW = +# FILE_PATTERNS = +# FORCE_LOCAL_INCLUDES = +# FORMULA_FONTSIZE = +# FORMULA_TRANSPARENT = +# FULL_PATH_NAMES = +# GENERATE_BUGLIST = +# GENERATE_DEPRECATEDLIST = +# GENERATE_LEGEND = +# GENERATE_TESTLIST = +# GENERATE_TODOLIST = +# GENERATE_TREEVIEW = +# GRAPHICAL_HIERARCHY = +# GROUP_GRAPHS = +# GROUP_NESTED_COMPOUNDS = +# HIDE_COMPOUND_REFERENCE= = +# HIDE_FRIEND_COMPOUNDS = +# HIDE_IN_BODY_DOCS = +# HIDE_SCOPE_NAMES = +# HIDE_UNDOC_CLASSES = +# HIDE_UNDOC_MEMBERS = +# HIDE_UNDOC_RELATIONS = +# HTML_COLORSTYLE_GAMMA = +# HTML_COLORSTYLE_HUE = +# HTML_COLORSTYLE_SAT = +# HTML_DYNAMIC_SECTIONS = +# HTML_EXTRA_FILES = +# HTML_EXTRA_STYLESHEET = +# HTML_FOOTER = +# HTML_HEADER = +# HTML_INDEX_NUM_ENTRIES = +# HTML_STYLESHEET = +# HTML_TIMESTAMP = +# IDL_PROPERTY_SUPPORT = +# IGNORE_PREFIX = +# IMAGE_PATH = +# INCLUDED_BY_GRAPH = +# INCLUDE_FILE_PATTERNS = +# INCLUDE_GRAPH = +# INCLUDE_PATH = +# INHERIT_DOCS = +# INLINE_GROUPED_CLASSES = +# INLINE_INFO = +# INLINE_INHERITED_MEMB = +# INLINE_SIMPLE_STRUCTS = +# INLINE_SOURCES = +# INPUT = +# INPUT_ENCODING = +# INTERACTIVE_SVG = +# INTERNAL_DOCS = +# JAVADOC_AUTOBRIEF = +# LAYOUT_FILE = +# MACRO_EXPANSION = +# MARKDOWN_SUPPORT = +# MAX_DOT_GRAPH_DEPTH = +# MSCFILE_DIRS = +# MULTILINE_CPP_IS_BRIEF = +# OPTIMIZE_FOR_FORTRAN = +# OPTIMIZE_OUTPUT_FOR_C = +# OPTIMIZE_OUTPUT_JAVA = +# OPTIMIZE_OUTPUT_VHDL = +# OUTPUT_LANGUAGE = +# PLANTUML_JAR_PATH = +# PREDEFINED = +# PROJECT_BRIEF = +# PROJECT_LOGO = +# PROJECT_NAME = +# PROJECT_NUMBER = +# QT_AUTOBRIEF = +# RECURSIVE = +# REFERENCED_BY_RELATION = +# REFERENCES_LINK_SOURCE = +# REFERENCES_RELATION = +# REPEAT_BRIEF = +# SEARCHENGINE = +# SEARCH_INCLUDES = +# SEPARATE_MEMBER_PAGES = +# SHORT_NAMES = +# SHOW_FILES = +# SHOW_GROUPED_MEMB_INC = +# SHOW_INCLUDE_FILES = +# SHOW_NAMESPACES = +# SHOW_USED_FILES = +# SIP_SUPPORT = +# SKIP_FUNCTION_MACROS = +# SORT_BRIEF_DOCS = +# SORT_BY_SCOPE_NAME = +# SORT_GROUP_NAMES = +# SORT_MEMBERS_CTORS_1ST = +# SORT_MEMBER_DOCS = +# SOURCE_BROWSER = +# SOURCE_TOOLTIPS = +# STRICT_PROTO_MATCHING = +# STRIP_CODE_COMMENTS = +# STRIP_FROM_INC_PATH = +# STRIP_FROM_PATH = +# SUBGROUPING = +# TAB_SIZE = +# TEMPLATE_RELATIONS = +# TREEVIEW_WIDTH = +# TYPEDEF_HIDES_STRUCT = +# UML_LIMIT_NUM_FIELDS = +# UML_LOOK = +# USE_MDFILE_AS_MAINPAGE = +# VERBATIM_HEADERS = +# \ No newline at end of file diff --git a/Doxyfile b/Doxyfile index a5fb8b3..bfd88ff 100644 --- a/Doxyfile +++ b/Doxyfile @@ -38,7 +38,7 @@ PROJECT_NAME = "DACE 2.0 API Manual" # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = "ESA 4000117860/16/NL/LF/as" +#PROJECT_NUMBER = "ESA 4000117860/16/NL/LF/as" # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/README.md b/README.md index b933d08..608e034 100644 --- a/README.md +++ b/README.md @@ -6,43 +6,32 @@ The Differential Algebra Computational Toolbox. ## Requirements To use the DACE library, you need a C++11 compatible C++ compiler. Almost all modern C++ compilers (GNU, Clang, MS Visual Studio) meet that requirement nowadays. We also highly recommend using CMake as a build system, and obviously Git as your source code management system. -On Windows, you need to download Microsoft Visual Studio (the free Community Edition is fine). On MacOS X, you need to install the Xcode package from Apple from the App Store. On Linux, you can install the compiler toolchain through your package system if they don't already come pre-installed. +On Linux, you can install the compiler toolchain through your package system if they don't already come pre-installed. On Windows it is suggested to use the WSL environment and follow the linux procedure, or download Microsoft Visual Studio (the free Community Edition is fine). On MacOS X, you need to install the Xcode package from Apple from the App Store. -To build Windows DACE library installer packages, you additionally need the NSIS installer compiler. +To build Windows DACE library installer packages, you additionally need the NSIS installer compiler. (Deprecated, will be removed soon) ## Getting started -To get started, either download one of the pre-built release packages for your platform (see the Releases tab above), or just clone the entire repository and build the library from scratch (which really is very easy). - -To use the DACE, we recommend using CMake. New DACE releases past version 2.0.1 come with a CMake package system plugin that allows very simple use of the DACE within the CMake framework. To get started, have a look at the Tutorials directories. Both of the Tutorials there are independent projects with their own stand-alone CMake file you can use as a starting point for your own programs. +To get started just clone the entire repository and build the library from scratch (which really is very easy). You can use the current development version if you want the latests improvement and bug fixes or point to the Releases page for official released versions. ## Building the DACE -To build your own version of the DACE library, simply clone this repository: +To build the DACE library, simply clone this repository: ``` -git clone "https://github.com/dacelib/dace.git" +git clone "https://github.com/dacelib/dace.git" dace ``` -Then create a build directory and run cmake, and then make to compile everything: +Then create a build directory, run cmake, and then cmake --build to compile everything: ``` -cd dace -mkdir build -cd build -cmake .. -make +mkdir dace-build +cmake -S dace/ -B dace-build/ +cmake --build dace-build/ ``` -After some compiling, you should have a sparkling new dace library ready for use. To install it directly into your system, you can just use: +After some compiling, you should have a sparkling new dace library ready for use in the dace-build folder. + +Optionally you can install the freshly built DACE library directly into your system: ``` -sudo make install +sudo cmake --install dace-build/ ``` The ```sudo``` is there to give you the required permissions to install into your system directories (usually ```/usr/local```). -## Running the Tutorials -To test if your installation was successful, you can try to build one of the Tutorials: -``` -cd ../Tutorials/Tutorial1 -mkdir build -cd build -cmake .. -make -``` -This should automatically compile all tutorials, ready to be run by you using e.g. the command ```./Example1```. +We have moved the tutorials to a different repository https://github.com/dacelib/dace-tutorials to keep the main repository clean. You can clone the tutorials repository and follow the instructions there to get started with DACE. -Also have a look at the [DACE Wiki pages](https://github.com/dacelib/dace/wiki) if you have further questions +Also have a look at the [DACE Wiki pages](https://github.com/dacelib/dace/wiki) if you have further questions. diff --git a/Tutorials/Tutorial1/CMakeLists.txt b/Tutorials/Tutorial1/CMakeLists.txt deleted file mode 100644 index f620300..0000000 --- a/Tutorials/Tutorial1/CMakeLists.txt +++ /dev/null @@ -1,67 +0,0 @@ -# DACE Examples -# -# This file shows how to use the DACE in a cmake project -# Simply import the DACE cmake package. It automatically provides the correct -# locations for headers and libraries. -# If the DACE is installed in a non-default location, tell cmake where to look: -# set(dace_DIR "/my/custom/root/usr/local/lib/cmake/dace") -# It exposes two targets: -# dace::dace (the dynamic DACE library) -# dace::dace_s (the static DACE library) -# Make sure to link your executable with one of those and include the correct -# header ( or ) in the code. -# -# On Windows only, the dace.dll dynamics library is typically expected to be -# located in the same directory as the executable using it. You can use a -# command such as the file(COPY ...) below to automatically copy the dace.dll -# from the DACE package next to your executables. - -cmake_minimum_required (VERSION 2.8.4) - -project(Examples1 CXX) - -find_package(dace 2.0.0 REQUIRED) - -if(WIN32) - get_target_property(DACEDLL dace::dace LOCATION) - file(COPY "${DACEDLL}" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}") -endif(WIN32) - -add_executable(Example1 Example1.cpp) -target_link_libraries(Example1 PUBLIC dace::dace) - -add_executable(Example2 Example2.cpp) -target_link_libraries(Example2 PUBLIC dace::dace) - -add_executable(Example3 Example3.cpp) -target_link_libraries(Example3 PUBLIC dace::dace) - -add_executable(Example4 Example4.cpp) -target_link_libraries(Example4 PUBLIC dace::dace) - -add_executable(Example5 Example5.cpp) -target_link_libraries(Example5 PUBLIC dace::dace) - -add_executable(Example6 Example6.cpp) -target_link_libraries(Example6 PUBLIC dace::dace) - -add_executable(Example7 Example7.cpp) -target_link_libraries(Example7 PUBLIC dace::dace) - -add_executable(Example8 Example8.cpp) -target_link_libraries(Example8 PUBLIC dace::dace) - -add_executable(Example9 Example9.cpp) -target_link_libraries(Example9 PUBLIC dace::dace) - -add_executable(Example10 Example10.cpp) -target_link_libraries(Example10 PUBLIC dace::dace) - -add_executable(Example11Pointwise Example11Pointwise.cpp) -target_link_libraries(Example11Pointwise PUBLIC dace::dace) - -add_executable(Example11DA Example11DA.cpp) -target_link_libraries(Example11DA PUBLIC dace::dace) - -add_executable(Example12 Example12.cpp) -target_link_libraries(Example12 PUBLIC dace::dace) diff --git a/Tutorials/Tutorial1/Example1.cpp b/Tutorials/Tutorial1/Example1.cpp deleted file mode 100644 index 2ea2c49..0000000 --- a/Tutorials/Tutorial1/Example1.cpp +++ /dev/null @@ -1,23 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - -int main( void ) -{ - - DA::init( 20, 1 ); // initialize DACE for 20th-order computations in 1 variable - - DA x = DA(1); // initialize x as DA - - DA y = sin(x); // compute y = sin(x) - - // print x and y to screen - cout << "x" << endl << x << endl; - cout << "y = sin(x)" << endl << y; - -} \ No newline at end of file diff --git a/Tutorials/Tutorial1/Example10.cpp b/Tutorials/Tutorial1/Example10.cpp deleted file mode 100644 index ac91a58..0000000 --- a/Tutorials/Tutorial1/Example10.cpp +++ /dev/null @@ -1,65 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -int main( void ) -{ - - DA::init( 10, 3 ); - - const double pi = 4.0*atan(1.0); - - // initialize cyl - AlgebraicVector cyl(3); - - cyl[0] = 100.0+DA(1); - cyl[1] = (0.0+DA(2))*pi/180.0; - cyl[2] = 0.0 + DA(3); - - // initialize cart and compute transformation - AlgebraicVector cart(3); - - cart[0] = cyl[0]*cos(cyl[1]); - cart[1] = cyl[0]*sin(cyl[1]);; - cart[2] = cyl[2]; - - // subtract constant part to build DirMap - AlgebraicVector DirMap(3); - - DirMap = cart-cart.cons(); - - cout << "Direct map: from cylindric to cartesian (DirMap)" << endl << endl; - cout << DirMap << endl << endl; - - // Invert DirMap to obtain InvMap - AlgebraicVector InvMap(3); - - InvMap = DirMap.invert(); - - cout << "Invers map: from cartesian to cylindric (InvMap)" << endl << endl; - cout << InvMap << endl << endl; - - getchar(); - - // Verification - cout << "Concatenate the direct and inverse map: (DirMap) o (InvMap) = DirMap(InvMap) = I" << endl << endl; - cout << DirMap.eval(InvMap) << endl; - -} - - - - - - - - - - - diff --git a/Tutorials/Tutorial1/Example11DA.cpp b/Tutorials/Tutorial1/Example11DA.cpp deleted file mode 100644 index a41c88f..0000000 --- a/Tutorials/Tutorial1/Example11DA.cpp +++ /dev/null @@ -1,104 +0,0 @@ -#include -#include -using namespace std; -using namespace DACE; - - -template AlgebraicVector TBP( AlgebraicVector x, double t ) -{ - - AlgebraicVector pos(3), res(6); - - pos[0] = x[0]; pos[1] = x[1]; pos[2] = x[2]; - - T r = pos.vnorm(); - - const double mu = 398600; // km^3/s^2 - - res[0] = x[3]; - res[1] = x[4]; - res[2] = x[5]; - - res[3] = -mu*pos[0]/(r*r*r); - res[4] = -mu*pos[1]/(r*r*r); - res[5] = -mu*pos[2]/(r*r*r); - - return res; - -} - - -template AlgebraicVector euler( AlgebraicVector x, double t0, double t1 ) -{ - const double hmax = 0.1; - int steps = ceil((t1-t0)/hmax); - double h = (t1-t0)/steps; - double t = t0; - - for( int i = 0; i < steps; i++ ) - { - x = x + h*TBP(x,t); - t += h; - } - - return x; -} - - -int main( void ) -{ - - DA::init( 3, 6 ); - - AlgebraicVector x0(6), xf(6); - - // Set initial conditions - const double mu = 398600; - const double ecc = 0.5; - - x0[0] = 6678.0 + DA(1); // 300 km altitude - x0[1] = 0.0 + DA(2); - x0[2] = 0.0 + DA(3); - x0[3] = 0.0 + DA(4); - x0[4] = sqrt(mu/6678.0)*sqrt(1+ecc) + DA(5); - x0[5] = 0.0 + DA(6); - - // integrate for half the orbital period - - const double pi = 4.0*atan(1.0); - - double a = 6678.0/(1-ecc); - - xf = euler( x0, 0.0, pi*sqrt(a*a*a/mu)); - - // print initial and final conditions - - cout << endl << "Initial conditions:" << endl << endl; - cout << x0 << endl << endl; - - cout << endl << "Final conditions:" << endl << endl; - cout << xf << endl << endl; - - getchar(); - - // Evaluate for a displaced initial condition - - AlgebraicVector Deltax0(6); - - Deltax0[0] = 1.0; // km - Deltax0[1] = -1.0; - Deltax0[2] = 0.0; - Deltax0[3] = 0.0; - Deltax0[4] = 0.0; - Deltax0[5] = 0.0; - - cout << "Displaced Initial condition:" << endl << endl; - cout << x0.cons() + Deltax0 << endl << endl; - - cout << "Displaced Final condition:" << endl << endl; - cout << xf.eval(Deltax0) << endl; - -} - - - diff --git a/Tutorials/Tutorial1/Example11DARK78.cpp b/Tutorials/Tutorial1/Example11DARK78.cpp deleted file mode 100644 index 340be00..0000000 --- a/Tutorials/Tutorial1/Example11DARK78.cpp +++ /dev/null @@ -1,256 +0,0 @@ -#include -#include -using namespace std; -using namespace DACE; - -#warning This example only works if the DACE was built with support for the (non-default) AlgebraicMatrix type! - -double min( double a, double b ) { - - double res = a; - if(a>b){res = b;} - return res; -} - - -double max( double a, double b ) { - - double res = a; - if(b>a){res = b;} - return res; -} - - -double normtmp( int N, vector X ) { - - int I; - double res = 0.0; - - for (I=0; I AlgebraicVector RK78( int N, AlgebraicVector Y0, double X0, double X1, AlgebraicVector (*f)(AlgebraicVector,double) ) -{ - - double ERREST; - double H0 = 0.001; double HS = 0.1; double H1 = 100.0; - double EPS = 1.e-12; double BS = 20*EPS; - - AlgebraicMatrix Z(N,16); - AlgebraicVector Y1(N); - vector Y1cons(N); - - double VIHMAX = 0.0, X, H; - int I, J, K; - double RFNORM, HH0, HH1; - - double HSQR = 1.0/9.0; - vector A(13), C(13), D(13); - AlgebraicMatrix B(13,12); - - A[0] = 0.0; A[1] = 1.0/18.0; A[2] = 1.0/12.0; A[3] = 1.0/8.0; A[4] = 5.0/16.0; A[5] = 3.0/8.0; - A[6] = 59.0/400.0; A[7] = 93.0/200.0; A[8] = 5490023248.0/9719169821.0; A[9] = 13.0/20.0; A[10] = 1201146811.0/1299019798.0; A[11] = 1.0; - A[12] = 1.0; - - B.at(0,0) = 0.0; B.at(0,1) = 0.0; B.at(0,2) = 0.0; B.at(0,3) = 0.0; B.at(0,4) = 0.0; - B.at(0,5) = 0.0; B.at(0,6) = 0.0; B.at(0,7) = 0.0; B.at(0,8) = 0.0; B.at(0,9) = 0.0; - B.at(0,10) = 0.0; B.at(0,11) = 0.0; - - B.at(1,0) = 1.0/18.0; B.at(1,1) = 0.0; B.at(1,2) = 0.0; B.at(1,3) = 0.0; B.at(1,4) = 0.0; - B.at(1,5) = 0.0; B.at(1,6) = 0.0; B.at(1,7) = 0.0; B.at(1,8) = 0.0; B.at(1,9) = 0.0; - B.at(1,10) = 0.0; B.at(1,11) = 0.0; - - B.at(2,0) = 1.0/48.0; B.at(2,1) = 1.0/16.0; B.at(2,2) = 0.0; B.at(2,3) = 0.0; B.at(2,4) = 0.0; - B.at(2,5) = 0.0; B.at(2,6) = 0.0; B.at(2,7) = 0.0; B.at(2,8) = 0.0; B.at(2,9) = 0.0; - B.at(2,10) = 0.0; B.at(2,11) = 0.0; - - B.at(3,0) = 1.0/32.0; B.at(3,1) = 0.0; B.at(3,2) = 3.0/32.0; B.at(3,3) = 0.0; B.at(3,4) = 0.0; - B.at(3,5) = 0.0; B.at(3,6) = 0.0; B.at(3,7) = 0.0; B.at(3,8) = 0.0; B.at(3,9) = 0.0; - B.at(3,10) = 0.0; B.at(3,11) = 0.0; - - B.at(4,0) = 5.0/16.0; B.at(4,1) = 0.0; B.at(4,2) = -75.0/64.0; B.at(4,3) = 75.0/64.0; B.at(4,4) = 0.0; - B.at(4,5) = 0.0; B.at(4,6) = 0.0; B.at(4,7) = 0.0; B.at(4,8) = 0.0; B.at(4,9) = 0.0; - B.at(4,10) = 0.0; B.at(4,11) = 0.0; - - B.at(5,0) = 3.0/80.0; B.at(5,1) = 0.0; B.at(5,2) = 0.0; B.at(5,3) = 3.0/16.0; B.at(5,4) = 3.0/20.0; - B.at(5,5) = 0.0; B.at(5,6) = 0.0; B.at(5,7) = 0.0; B.at(5,8) = 0.0; B.at(5,9) = 0.0; - B.at(5,10) = 0.0; B.at(5,11) = 0.0; - - B.at(6,0) = 29443841.0/614563906.0; B.at(6,1) = 0.0; B.at(6,2) = 0.0; B.at(6,3) = 77736538.0/692538347.0; B.at(6,4) = -28693883.0/1125000000.0; - B.at(6,5) = 23124283.0/1800000000.0; B.at(6,6) = 0.0; B.at(6,7) = 0.0; B.at(6,8) = 0.0; B.at(6,9) = 0.0; - B.at(6,10) = 0.0; B.at(6,11) = 0.0; - - B.at(7,0) = 16016141.0/946692911.0; B.at(7,1) = 0.0; B.at(7,2) = 0.0; B.at(7,3) = 61564180.0/158732637.0; B.at(7,4) = 22789713.0/633445777.0; - B.at(7,5) = 545815736.0/2771057229.0; B.at(7,6) = -180193667.0/1043307555.0; B.at(7,7) = 0.0; B.at(7,8) = 0.0; B.at(7,9) = 0.0; - B.at(7,10) = 0.0; B.at(7,11) = 0.0; - - B.at(8,0) = 39632708.0/573591083.0; B.at(8,1) = 0.0; B.at(8,2) = 0.0; B.at(8,3) = -433636366.0/683701615.0; B.at(8,4) = -421739975.0/2616292301.0; - B.at(8,5) = 100302831.0/723423059.0; B.at(8,6) = 790204164.0/839813087.0; B.at(8,7) = 800635310.0/3783071287.0; B.at(8,8) = 0.0; B.at(8,9) = 0.0; - B.at(8,10) = 0.0; B.at(8,11) = 0.0; - - B.at(9,0) = 246121993.0/1340847787.0; B.at(9,1) = 0.0; B.at(9,2) = 0.0; B.at(9,3) = -37695042795.0/15268766246.0; B.at(9,4) = -309121744.0/1061227803.0; - B.at(9,5) = -12992083.0/490766935.0; B.at(9,6) = 6005943493.0/2108947869.0; B.at(9,7) = 393006217.0/1396673457.0; B.at(9,8) = 123872331.0/1001029789.0; B.at(9,9) = 0.0; - B.at(9,10) = 0.0; B.at(9,11) = 0.0; - - B.at(10,0) = -1028468189.0/846180014.0; B.at(10,1) = 0.0; B.at(10,2) = 0.0; B.at(10,3) = 8478235783.0/508512852.0; B.at(10,4) = 1311729495.0/1432422823.0; - B.at(10,5) = -10304129995.0/1701304382.0; B.at(10,6) = -48777925059.0/3047939560.0; B.at(10,7) = 15336726248.0/1032824649.0; B.at(10,8) = -45442868181.0/3398467696.0; B.at(10,9) = 3065993473.0/597172653.0; - B.at(10,10) = 0.0; B.at(10,11) = 0.0; - - B.at(11,0) = 185892177.0/718116043.0; B.at(11,1) = 0.0; B.at(11,2) = 0.0; B.at(11,3) = -3185094517.0/667107341.0; B.at(11,4) = -477755414.0/1098053517.0; - B.at(11,5) = -703635378.0/230739211.0; B.at(11,6) = 5731566787.0/1027545527.0; B.at(11,7) = 5232866602.0/850066563.0; B.at(11,8) = -4093664535.0/808688257.0; B.at(11,9) = 3962137247.0/1805957418.0; - B.at(11,10) = 65686358.0/487910083.0; B.at(11,11) = 0.0; - - B.at(12,0) = 403863854.0/491063109.0; B.at(12,1) = 0.0; B.at(12,2) = 0.0; B.at(12,3) = - 5068492393.0/434740067.0; B.at(12,4) = -411421997.0/543043805.0; - B.at(12,5) = 652783627.0/914296604.0; B.at(12,6) = 11173962825.0/925320556.0; B.at(12,7) = -13158990841.0/6184727034.0; B.at(12,8) = 3936647629.0/1978049680.0; B.at(12,9) = -160528059.0/685178525.0; - B.at(12,10) = 248638103.0/1413531060.0; B.at(12,11) = 0.0; - - C[0] = 14005451.0/335480064.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; C[4] = 0.0; C[5] = -59238493.0/1068277825.0; - C[6] = 181606767.0/758867731.0; C[7] = 561292985.0/797845732.0; C[8] = -1041891430.0/1371343529.0; C[9] = 760417239.0/1151165299.0; C[10] = 118820643.0/751138087.0; C[11] = -528747749.0/2220607170.0; - C[12] = 1.0/4.0; - - D[0] = 13451932.0/455176623.0; D[1] = 0.0; D[2] = 0.0; D[3] = 0.0; D[4] = 0.0; D[5] = -808719846.0/976000145.0; - D[6] = 1757004468.0/5645159321.0; D[7] = 656045339.0/265891186.0; D[8] = -3867574721.0/1518517206.0; D[9] = 465885868.0/322736535.0; D[10] = 53011238.0/667516719.0; D[11] = 2.0/45.0; - D[12] = 0.0; - - for( I = 0; I < N; I++ ) - { - Z.at(I,0) = Y0[I]; - Z.at(I,1) = 0.0 ; - } - - H = abs(HS) ; HH0 = abs(H0) ; HH1 = abs(H1) ; - X = X0 ; RFNORM = 0.0 ; ERREST = 0.0 ; - - while(X != X1){ - - // compute new stepsize - if (RFNORM != 0) {H = H*min(4.0,exp(HSQR*log(EPS/RFNORM)));} - if (abs(H)>abs(HH1)) { H = HH1; } else if ( abs(H)0) { H = X1-X; } - - for (J = 0; J<13; J++) { - - for (I = 0; IBS) && (abs(H/H0)>1.2)) { - H = H/3.0; - RFNORM = 0; - } - else { - for (I = 0; I AlgebraicVector rhs( AlgebraicVector x, double t ) -{ - - AlgebraicVector pos(3), res(6); - - pos[0] = x[0]; pos[1] = x[1]; pos[2] = x[2]; - - T r = pos.vnorm(); - - const double mu = 398600; // km^3/s^2 - - res[0] = x[3]; - res[1] = x[4]; - res[2] = x[5]; - - res[3] = -mu*pos[0]/(r*r*r); - res[4] = -mu*pos[1]/(r*r*r); - res[5] = -mu*pos[2]/(r*r*r); - - return res; - -} - - -int main( void ) -{ - - DA::init( 3, 6 ); - - AlgebraicVector x0(6), xf(6); - - // Set initial conditions - const double mu = 398600; - const double ecc = 0.5; - - x0[0] = 6678.0 + DA(1); // 300 km altitude - x0[1] = 0.0 + DA(2); - x0[2] = 0.0 + DA(3); - x0[3] = 0.0 + DA(4); - x0[4] = sqrt(mu/6678.0)*sqrt(1+ecc) + DA(5); - x0[5] = 0.0 + DA(6); - - // integrate for half the orbital period - - const double pi = 4.0*atan(1.0); - - double a = 6678.0/(1-ecc); - - xf = RK78( 6, x0, 0.0, pi*sqrt(a*a*a/mu), rhs ); - - // print initial and final conditions - - cout << endl << "Initial conditions:" << endl << endl; - cout << x0 << endl << endl; - - cout << endl << "Final conditions:" << endl << endl; - cout << xf << endl; - -} - - - diff --git a/Tutorials/Tutorial1/Example11Pointwise.cpp b/Tutorials/Tutorial1/Example11Pointwise.cpp deleted file mode 100644 index d0f1fc7..0000000 --- a/Tutorials/Tutorial1/Example11Pointwise.cpp +++ /dev/null @@ -1,85 +0,0 @@ -#include -#include -using namespace std; -using namespace DACE; - - -template AlgebraicVector TBP( AlgebraicVector x, double t ) -{ - - AlgebraicVector pos(3), res(6); - - pos[0] = x[0]; pos[1] = x[1]; pos[2] = x[2]; - - T r = pos.vnorm(); - - const double mu = 398600; // km^3/s^2 - - res[0] = x[3]; - res[1] = x[4]; - res[2] = x[5]; - - res[3] = -mu*pos[0]/(r*r*r); - res[4] = -mu*pos[1]/(r*r*r); - res[5] = -mu*pos[2]/(r*r*r); - - return res; - -} - - -template AlgebraicVector euler( AlgebraicVector x, double t0, double t1 ) -{ - const double hmax = 0.1; - int steps = ceil((t1-t0)/hmax); - double h = (t1-t0)/steps; - double t = t0; - - for( int i = 0; i < steps; i++ ) - { - x = x + h*TBP(x,t); - t += h; - } - - return x; -} - - - - -int main( void ) -{ - - AlgebraicVector x0(6), xf(6); - - // Set initial conditions - const double mu = 398600; - const double ecc = 0.5; - - x0[0] = 6678.0; // 300 km altitude - x0[1] = 0.0; - x0[2] = 0.0; - x0[3] = 0.0; - x0[4] = sqrt(mu/6678.0)*sqrt(1+ecc); - x0[5] = 0.0; - - // integrate for half the orbital period - - const double pi = 4.0*atan(1.0); - - double a = 6678.0/(1-ecc); - - xf = euler( x0, 0.0, pi*sqrt(a*a*a/mu)); - - // print initial and final conditions - - cout << endl << "Initial conditions:" << endl << endl; - cout << x0 << endl << endl; - - cout << endl << "Final conditions:" << endl << endl; - cout << xf << endl; - -} - - - diff --git a/Tutorials/Tutorial1/Example12.cpp b/Tutorials/Tutorial1/Example12.cpp deleted file mode 100644 index ed090ab..0000000 --- a/Tutorials/Tutorial1/Example12.cpp +++ /dev/null @@ -1,80 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -int main( void ) -{ - - DA::init( 10, 2 ); - - const double tol = 1.0e-12; - - const double pi = 4.0*atan(1.0); - const double mu = 1.0; - - DA a = 1.0; - DA e = 0.5; - double t = pi/2.0; - - DA M = sqrt(mu/(a*a*a))*t; // real at this stage - - DA EccAn = M; //first guess - - double err = abs(EccAn - e*sin(EccAn) - M); - - // Newton's method for the reference solution - while(err>tol) - { - EccAn = EccAn - (EccAn - e*sin(EccAn) - M)/(1 - e*cos(EccAn)); - - err = abs(EccAn - e*sin(EccAn) - M); - - } - - cout << setprecision(12); - cout << "Reference solution: E = " << cons(EccAn) << endl; - - getchar(); - - a = 1.0 + DA(1); - e = 0.5 + DA(2); - - M = sqrt(mu/(a*a*a))*t; // now M is a DA - - // Newton's method for the Taylor expansion of the solution - int i = 1; - while( i <= 10){ - - EccAn = EccAn - (EccAn - e*sin(EccAn) - M)/(1 - e*cos(EccAn)); - i *= 2; - - } - - cout << "Taylor expansion of E" << cons(EccAn) << endl; - cout << EccAn << endl; - - getchar(); - - cout << "Let's verify it is the Taylor expansion of the solution:" << endl; - cout << "Evaluate (E - e*sin(E) - M) in DA" << endl; - - cout << EccAn - e*sin(EccAn) - M << endl; - -} - - - - - - - - - - - diff --git a/Tutorials/Tutorial1/Example2.cpp b/Tutorials/Tutorial1/Example2.cpp deleted file mode 100644 index 7f84f89..0000000 --- a/Tutorials/Tutorial1/Example2.cpp +++ /dev/null @@ -1,28 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - -int main( void ) -{ - - DA::init( 20, 1 ); - - DA x = DA(1); - - DA y1 = sqr(sin(x)); // compute sin(x)^2 - - cout << "sin(x)^2" << endl << y1 << endl; - - DA y2 = sqr(cos(x)); // compute cos(x)^2 - - cout << "cos(x)^2" << endl << y2 << endl; - - // compute and print sin(x)^2+cos(x)^2 - cout << "sin(x)^2+cos(x)^2" << endl << y1+y2 << endl; - -} \ No newline at end of file diff --git a/Tutorials/Tutorial1/Example3.cpp b/Tutorials/Tutorial1/Example3.cpp deleted file mode 100644 index 5ea64cf..0000000 --- a/Tutorials/Tutorial1/Example3.cpp +++ /dev/null @@ -1,24 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -int main( void ) -{ - - DA::init( 1, 1 ); // Initialize DACE for 1st-order computations in 1 variable - - DA x = 3 + DA(1); // Initialize x as DA around 3 - - cout << "x" << endl << x << endl; - - DA f = 1/(x+1/x); // Evaluate f(x) = 1/(x+1/x) - - cout << "f(x) = 1/(x+1/x)" << endl << f << endl; - -} \ No newline at end of file diff --git a/Tutorials/Tutorial1/Example4.cpp b/Tutorials/Tutorial1/Example4.cpp deleted file mode 100644 index 8c18dc8..0000000 --- a/Tutorials/Tutorial1/Example4.cpp +++ /dev/null @@ -1,30 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -int main( void ) -{ - - DA::init( 20, 1 ); - - DA x = DA(1); - - DA y = cos(x)-1; // Compute [cos(x)-1] - - cout << "[cos(x)-1]" << endl << y << endl; - - // Compute [cos(x)-1]^11 - for ( int i = 0; i < 10; i++) - { - y = y*(cos(x)-1); - } - - cout << "[cos(x)-1]^11" << endl << y << endl; - -} \ No newline at end of file diff --git a/Tutorials/Tutorial1/Example5.cpp b/Tutorials/Tutorial1/Example5.cpp deleted file mode 100644 index 31d0748..0000000 --- a/Tutorials/Tutorial1/Example5.cpp +++ /dev/null @@ -1,34 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -int main( void ) -{ - - DA::init( 20, 1 ); - - DA x = DA(1); - - DA y = sin(x); - - DA dy = y.deriv(1); // compute Taylor expansion of d[sin(x)]/dx - - // print d[sin(x)]/dx and cos(x) to compare - cout << "d[sin(x)]/dx" << endl << dy << endl; - cout << "cos(x)" << endl << cos(x) << endl; - - getchar(); - - DA int_y = y.integ(1); // compute Taylor expansion of int[sin(x)dx] - - // print int[sin(x)dx] and -cos(x) to compare - cout << "int[sin(x)dx]" << endl << int_y << endl; - cout << "-cos(x)" << endl << -cos(x) << endl; - -} \ No newline at end of file diff --git a/Tutorials/Tutorial1/Example6.cpp b/Tutorials/Tutorial1/Example6.cpp deleted file mode 100644 index 3fd8be5..0000000 --- a/Tutorials/Tutorial1/Example6.cpp +++ /dev/null @@ -1,38 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -template T ErrFunc(T x) -{ - const double pi = 4.0*atan(1.0); - - T z = 1.0/sqrt(2.0*pi)*exp(-x*x/2); - return z; -} - - -int main( void ) -{ - - DA::init( 24, 1 ); // initialize DACE for 20th-order computations in 1 variable - - DA x = DA(1); - - DA y = ErrFunc(x); // compute Taylor expansion of the erf - - DA Inty = y.integ(1); // compute the Taylor expansion of the indefinite integral of erf - - double value = Inty.evalScalar(1.0)-Inty.evalScalar(-1.0); // compute int_{-1}^{+1} (erf) - - cout << "int_{-1}^{+1} (erf)" << endl; - cout << setprecision(12); - cout << "Exact result: " << 0.682689492137 << endl; - cout << "Approx result " << value << endl; - -} \ No newline at end of file diff --git a/Tutorials/Tutorial1/Example7.cpp b/Tutorials/Tutorial1/Example7.cpp deleted file mode 100644 index 3279c02..0000000 --- a/Tutorials/Tutorial1/Example7.cpp +++ /dev/null @@ -1,52 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -template T somb(AlgebraicVector x) -{ - - return sin(sqrt(sqr(x[0]) + sqr(x[1])))/sqrt(sqr(x[0]) + sqr(x[1])); - -} - - -int main( void ) -{ - - DA::init( 10, 2 ); // initialize DACE for 10th-order computations in 2 variables - - cout << "Initialize x as two-dim DA vector around (0,0)" << endl << endl; - - AlgebraicVector x(2); - - x[0] = DA(1); - x[1] = DA(2); - - cout << "x" << endl << x << endl; - - getchar(); - - cout << "Evaluate sombrero function" << endl << endl; - - DA z = somb(x); // Evaluate sombrero function - - cout << "z" << endl << z << endl; - -} - - - - - - - - - - - diff --git a/Tutorials/Tutorial1/Example8.cpp b/Tutorials/Tutorial1/Example8.cpp deleted file mode 100644 index 5da17f0..0000000 --- a/Tutorials/Tutorial1/Example8.cpp +++ /dev/null @@ -1,54 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -template T somb(AlgebraicVector x) -{ - - return sin(sqrt(sqr(x[0]) + sqr(x[1])))/sqrt(sqr(x[0]) + sqr(x[1])); - -} - - -int main( void ) -{ - - DA::init( 1, 2 ); // initialize DACE for 1st-order computations in 2 variables - - cout << "Initialize x as two-dim DA vector around (2,3)" << endl << endl; - - AlgebraicVector x(2); - - x[0] = 2.0 + DA(1); - x[1] = 3.0 + DA(2); - - DA z = somb(x); // Evaluate sombrero function - - cout << "Sombrero function" << endl << endl; - cout << z << endl << endl; - - AlgebraicVector grad_z(2); // Declare DA vector that will contain the gradient of sombrero function - - grad_z = z.gradient(); // Compute gradient of sombrero function - - cout << "Gradient of sombrero function" << endl << endl; - cout << grad_z << endl; - -} - - - - - - - - - - - diff --git a/Tutorials/Tutorial1/Example9.cpp b/Tutorials/Tutorial1/Example9.cpp deleted file mode 100644 index a3563ec..0000000 --- a/Tutorials/Tutorial1/Example9.cpp +++ /dev/null @@ -1,42 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace DACE; - - -int main( void ) -{ - - DA::init( 10, 1 ); - - DA x = DA(1); - - AlgebraicVector y(1), inv_y(1); - - y[0] = sin(x); // Compute Taylor expansion of sin(x) - - inv_y = y.invert(); // Invert Taylor polynomial - - // Compare with asin(x) - cout << "Polynomial inversion of sin(x)" << endl << endl; - cout << inv_y << endl << endl; - - cout << "asin(x)" << endl << endl; - cout << asin(x) << endl; - -} - - - - - - - - - - - diff --git a/Tutorials/Tutorial1/Tutorial1.pdf b/Tutorials/Tutorial1/Tutorial1.pdf deleted file mode 100644 index 69cc688..0000000 Binary files a/Tutorials/Tutorial1/Tutorial1.pdf and /dev/null differ diff --git a/Tutorials/Tutorial2/1Basics-Ex.cpp b/Tutorials/Tutorial2/1Basics-Ex.cpp deleted file mode 100644 index ec6a766..0000000 --- a/Tutorials/Tutorial2/1Basics-Ex.cpp +++ /dev/null @@ -1,74 +0,0 @@ -#include -#include - -using namespace std; -using namespace DACE; - -// Exercise 1.1.2: first steps -void ex1_1_2( ) -{ - DA x = DA(1); - DA func = 3*(x+3)-x+1-(x+8); - cout << "Exercise 1.1.2: First steps\n" << func << endl; -} - -// Exercise 1.1.3: different expansion point -void ex1_1_3( ) -{ - DA x = DA(1); - DA func = sin(1.0+x); - cout << "Exercise 1.1.3: Different expansion point\n" << func << endl; -} - -// Exercise 1.1.4: a higer power -void ex1_1_4( ) -{ - DA x = DA(1); - DA func = sin(x); - DA res = 1.0; // this makes res a constant function P(x) = 1.0 - for( int i = 0; i < 11; i++ ) res = res*sin(x); - cout << "Exercise 1.1.4: A higher power\n" << res << endl; -} - -// Exercise 1.1.5: two arguments -template T ex1_1_5( T x, T y ) -{ - return sqrt( 1.0 + x*x + y*y ); -} - -// Exercise 1.2.1: identity crisis -void ex1_2_1( ) -{ - DA x = DA(1); - DA s2 = sin(x)*sin(x); - DA c2 = cos(x)*cos(x); - cout << "Exercise 1.2.1: Identity crisis\n" << s2 << c2 << s2+c2 << endl; -} - -// Exercise 1.2.2: Breaking bad -template T ex1_2_2( T x, T y ) -{ - T r = sqrt( x*x + y*y ); - return sin(r)/r; -} - -int main( void ) -{ - DA::init( 10, 2 ); - cout.precision( 15 ); - - DA x = DA(1); - DA y = DA(2); - - ex1_1_2( ); - ex1_1_3( ); - ex1_1_4( ); - cout << "Exercise 1.1.5: Two arguments\n" << ex1_1_5( x, y ) << endl; - - ex1_2_1( ); - try { - cout << "Exercise 1.2.2: Breaking bad\n" - << ex1_2_2( 0.0, 0.0 ) << endl - << ex1_2_2( x, y ) << endl; - } catch( exception const& ex ) { cout << ex.what( ) << endl; } -} \ No newline at end of file diff --git a/Tutorials/Tutorial2/2IntDiff-Ex.cpp b/Tutorials/Tutorial2/2IntDiff-Ex.cpp deleted file mode 100644 index 3741669..0000000 --- a/Tutorials/Tutorial2/2IntDiff-Ex.cpp +++ /dev/null @@ -1,114 +0,0 @@ -#include -#include -using namespace std; using namespace DACE; - -const double pi = 3.141592653589793; - -// Exercise 2.1.1: derivatives -template T somb( T x, T y ) -{ -// return sin(y)*exp(x*x); - T r = sqrt( x*x + y*y ); - return sin(r)/r; -} - -void ex2_1_1( ) -{ - double x0 = 2.0, y0 = 3.0; - DA x = DA(1); DA y = DA(2); - DA func = somb( x0+x, y0+y ); // expand sombrero function around (x0,y0) - - // compute the derivative using DA - double dadx = cons(func.deriv(1)); - double dady = cons(func.deriv(2)); - double dadxx = cons(func.deriv(1).deriv(1)); - double dadxy = cons(func.deriv(1).deriv(2)); - double dadyy = cons(func.deriv(2).deriv(2)); - double dadxxx = cons(func.deriv(1).deriv(1).deriv(1)); - - // compute the derivatives using divided differences - const double h = 1e-3; - double dx = (somb(x0+h,y0)-somb(x0-h,y0))/(2.0*h); - double dy = (somb(x0,y0+h)-somb(x0,y0-h))/(2.0*h); - double dxx = (somb(x0+2.0*h,y0)-2.0*somb(x0,y0)+somb(x0-2.0*h,y0))/(4.0*h*h); - double dxy = (somb(x0+h,y0+h)-somb(x0-h,y0+h)-somb(x0+h,y0-h)+somb(x0-h,y0-h))/(4.0*h*h); - double dyy = (somb(x0,y0+2.0*h)-2.0*somb(x0,y0)+somb(x0,y0-2.0*h))/(4.0*h*h); - double dxxx = (somb(x0+3.0*h,y0)-3.0*somb(x0+h,y0)+3.0*somb(x0-h,y0)-somb(x0-3.0*h,y0))/(8.0*h*h*h); - - cout << "Exercise 2.1.1: Numerical derivatives\n" - << "d/dx: " << abs(dadx-dx) << endl - << "d/dy: " << abs(dady-dy) << endl - << "d/dxx: " << abs(dadxx-dxx) << endl - << "d/dxy: " << abs(dadxy-dxy) << endl - << "d/dyy: " << abs(dadyy-dyy) << endl - << "d/dxxx: " << abs(dadxxx-dxxx) << endl << endl; -} - -// Exercise 2.1.2: indefinite integral -void ex2_1_2( ) -{ - DA x = DA(1); - DA func = (1.0/(1+sqr(x))).integ(1); // DA integral - DA integral = atan(x); // analytical integral DA expanded - cout << "Exercise 2.1.2: Indefinite integral\n" << func-integral << endl; -} - -// Exercise 2.1.3: expand the error function -void ex2_1_3( ) -{ - DA t = DA(1); - DA erf = 2.0/sqrt(pi)*exp(-sqr(t)).integ(1); // error function erf(x) - cout << "Exercise 2.1.3: Error function\n" << erf << endl; -} - -// Exercise 2.2.1: DA based Newton solver -template T f( T x ) -{ - return x*sin(x)+cos(x); -// const double e = 0.1; const double M = 1.0; -// return x-e*sin(x)-M; // Kepler's equation -} - -double ex2_2_1( double x0 ) -{ - DA::pushTO( 1 ); // for this Newton solver we only need first derivatives - - const double err = 1e-14; - DA x = DA(1), func; - int i = 0; - - do { - func = f( x0+x ); // expand f around x0 - x0 = x0 - cons(func)/cons(func.deriv(1)); // Newton step - i++; - } while( (abs(cons(func))>err) && (i<1000) ); - - cout << "Exercise 2.2.1: DA Newton solver\n" - << "root x0: " << x0 << endl - << "residue at f(x0): " << abs(f(x0)) << endl - << "Newton iterations: " << i << endl << endl; - - DA::popTO( ); // don't forget to reset computation order to old value for following computations - return x0; -} - -// Exercise 2.2.2: expand the error function around x0 -void ex2_2_2( double x0 ) -{ - DA t = x0+DA(1); - DA erf = 2.0/sqrt(pi)*exp(-sqr(t)).integ(1); // error function erf(x) - cout << "Exercise 2.2.2: Shifted indefinite integral\n" << erf << endl; -} - -int main( void ) -{ - DA::init( 30, 2 ); - cout.precision( 15 ); - - ex2_1_1( ); - ex2_1_2( ); - ex2_1_3( ); - - ex2_2_1( 3.6 ); - ex2_2_2( 1.0 ); -} \ No newline at end of file diff --git a/Tutorials/Tutorial2/3Evaluation-Ex.cpp b/Tutorials/Tutorial2/3Evaluation-Ex.cpp deleted file mode 100644 index 66b05af..0000000 --- a/Tutorials/Tutorial2/3Evaluation-Ex.cpp +++ /dev/null @@ -1,154 +0,0 @@ -#include -#include -#include -using namespace std; using namespace DACE; - -const double pi = 3.141592653589793; -const double val1 = 1.685401585899429; // reference value of the integral over -1, +1 -const double val2 = 1.990644530037905; // reference value of the integral over -2, +2 - -// Exercise 3.1.1: plot a 1D polynomial -void ex3_1_1( ) -{ - const double x0 = 0.0; // expansion point - const int N = 100; // number of points in grid - const double hw = 2.0; // length of grid in each direction from expansion point - - DA x = DA(1)+x0; - DA func = exp(-x*x); - ofstream file; - - file.open( "ex3_1_1.dat" ); - for( int i = 0; i < N; i++ ) - { - double xx = -hw + i*2.0*hw/(N-1); // point on the grid on [-hw,hw] - double rda = func.evalScalar( xx ); // note: this is not an efficient way to repeatedly evaluate the same polynomial - xx += x0; // add expansion point x0 for double evaluation - double rdouble = exp(-xx*xx); - - file << xx << " " << rda << " " << rdouble << endl; - } - file.close( ); - // gnuplot command: plot 'ex3_1_1.dat'u 1:2 w l, 'ex3_1_1.dat' u 1:3 w l - // or for the error: plot 'ex3_1_1.dat' u 1:($2-$3) w l -} - -// Exercise 3.1.2: plot a 2D polynomial -template T somb( T x, T y ) -{ - T r = sqrt( x*x + y*y ); - return sin(r)/r; -} - -void ex3_1_2( ) -{ - const double x0 = 1.0; // expansion point x - const double y0 = 1.0; // expansion point y - const int N = 50; // number of points in grid - - DA x = DA(1)+x0; DA y = DA(2)+y0; - DA func = somb(x,y); - vector arg(2); // vector holding two doubles - ofstream file; - - file.open( "ex3_1_2.dat" ); - for( int i = 0; i < N; i++ ) - { - arg[0] = -1.0 + i*2.0/(N-1); // x coordinate on the grid on [-1,1] - for( int j = 0; j < N; j++ ) - { - arg[1] = -1.0 + j*2.0/(N-1); // y coordinate on the grid on [-1,1] - double rda = func.eval( arg ); // note: this is not an efficient way to repeatedly evaluate the same polynomial - double rdouble = somb( x0+arg[0], y0+arg[1] ); - file << arg[0] << " " << arg[1] << " " << rda << " " << rdouble << endl; - } - file << endl; // empty line between lines of data for gnuplot - } - file.close( ); - // gnuplot command: splot 'ex3_1_2.dat' u 1:2:3 w l, 'ex3_1_2.dat' u 1:2:4 w l - // or for the error: splot 'ex3_1_2.dat' u 1:2:($3-$4) w l -} - -// Exercise 3.1.3: Sinusitis -void ex3_1_3( ) -{ - DA::pushTO( 10 ); - - DA x = DA(1); - DA sinda = sin(x); - DA res1 = sin( x+2 ); // compute directly sin(1+DA(1)) - DA res2 = sinda.evalScalar( x+2 ); // evaluate expansion of sine with 1+DA(1) - - cout << "Exercise 4.1.3: Sinusitis" << endl << res1-res2 << endl; - - DA::popTO( ); -} - -// Exercise 3.1.4: Gauss integral I -void ex3_1_4( void ) -{ - ofstream file; - - file.open( "ex3_1_4.dat" ); - file.precision( 16 ); // print 16 digits, similar to MATLAB's "format long" - for( int order = 1; order <= 40; order++ ) - { - DA::setTO( order ); // limit the computation order - DA t = DA(1); - DA erf = 2.0/sqrt(pi)*exp(-sqr(t)).integ(1); // error function erf(x) - double res = erf.evalScalar( 1.0 )-erf.evalScalar( -1.0 ); - file << order << " " << res << " " << log(abs(res-val1))/log(10.0) << endl; - } - file.close( ); - // gnuplot command: plot 'ex3_1_4.dat'u 1:2 w l - // or for the error: plot 'ex3_1_4.dat'u 1:3 w l -} - -// Exercise 3.2.1 & 3.2.2: Gauss integral II -double gaussInt( double a, double b ) // compute integral of Gaussian on interval [a,b] -{ - DA t = (a+b)/2.0+DA(1); // expand around center point - DA func = 2.0/sqrt(pi)*exp(-sqr(t)).integ(1); - return func.evalScalar( (b-a)/2.0 ) - func.evalScalar( -(b-a)/2.0 ); // evaluate over -+ half width -} - -void ex3_2_1( void ) -{ - const double hw = 2.0; // half-width of the interval to integrate on, i.e. [-hw,hw] - ofstream file; - double res; - - file.open( "ex3_2_1.dat" ); - file.precision( 16 ); - DA::pushTO( 9 ); - for( int N = 1; N <= 30; N++ ) - { - res = 0.0; - for( int i = 1; i <= N; i++ ) - { - double ai = -hw + (i-1)*2.0*hw/N; - double ai1 = -hw + i*2.0*hw/N; - res += gaussInt( ai, ai1 ); - } - file << N << " " << res << " " << log(abs(res-val2))/log(10.0) << endl; - } - DA::popTO( ); - // compare to single expansion at full computation order - res = gaussInt( -hw, hw ); - file << endl << 1 << " " << res << " " << log(abs(res-val2))/log(10.0) << endl; - file.close( ); - // gnuplot command: plot 'ex3_2_1.dat'u 1:3 w lp -} - -int main( void ) -{ - DA::init( 30, 2 ); // init with maximum computation order - - ex3_1_1( ); - ex3_1_2( ); - ex3_1_3( ); - ex3_1_4( ); - - ex3_2_1( ); -} - diff --git a/Tutorials/Tutorial2/4Solver-Ex.cpp b/Tutorials/Tutorial2/4Solver-Ex.cpp deleted file mode 100644 index cfae03b..0000000 --- a/Tutorials/Tutorial2/4Solver-Ex.cpp +++ /dev/null @@ -1,138 +0,0 @@ -#include -#include -#include -using namespace std; using namespace DACE; - -const int order = 10; - -template T Nf( T x, T p ) -{ - return x - (x*x-p)/(2*x); -} - -// Exercise 4.1.3: Naive Newton -void ex4_1_3( ) -{ - double tol = 1e-14; // tolerance - double p0 = 4; double x0 = 1; // expansion point and initial guess - DA p = p0 + DA(1); // DA parameter - DA x = x0; // DA initial guess - DA xp; - int i = 0; - - do - { - xp = x; - x = Nf( xp, p ); - i++; - } while( (abs(xp-x) > tol) && (i<1000) ); - cout << "Exercise 4.1.3: Naive Newton" << endl << x << endl << sqrt(p)-x - << "Number of iterations: " << i << endl << endl; -} - -// Exercise 4.1.4: complicated parameters -void ex4_1_4( ) -{ - double p0 = 0; double x0 = 1; // x0 must now satisfy f(x0,cos(p0))=0 - DA p = cos( p0+DA(1) ); - DA x = x0; - - int i = 1; - while( i <= order ) - { - x = Nf( x, p ); - i *= 2; - } - cout << "Exercise 4.1.3: Naive Newton" << endl << x << endl << sqrt(p)-x; -} - -// Exercise 4.2.1: Full DA Newton solver -void ex4_2_1( double p0 ) -{ - const double tol = 1e-14; - double x0 = p0/2.0, xp; // x0 is just some initial guess - int i = 0; - - // double precision computation => fast - do{ - xp = x0; - x0 = Nf( xp, p0 ); - i++; - } while( (abs(xp-x0) > tol) && (i<1000) ); - - // DA computation => slow - DA p = p0 + DA(1); - DA x = x0; - i = 1; - while( i <= order ) - { - x = Nf( x, p ); - i *= 2; - } - - cout << "Exercise 4.2.1: Full DA Newton" << endl << x << endl << sqrt(p)-x; -} - -// Exercise 4.2.2 & 4.2.3: Kepler's equation solver -template T NKep( T E, T M, T ecc ) -{ - return E - (E-ecc*sin(E)-M)/(1-ecc*cos(E)); -} - -// double precision Kepler solver -double Kepler( double M, double ecc ) -{ - const double tol = 1e-14; - double E0 = M, Ep; - int i = 0; - - do{ - Ep = E0; - E0 = NKep( Ep, M, ecc ); - i++; - } while( (abs(Ep-E0) > tol) && (i<1000) ); - - return E0; -} - -void ex4_2_2( double M0, double ecc0 ) -{ - const double tol = 1e-14; - - DA M = M0 + DA(1); - DA E = Kepler( M0, ecc0 ); // reference solution - DA ecc = ecc0; // keep eccentricity constant (4.2.2) - //DA ecc = ecc0 + 0.1*DA(2); // also expand w.r.t. eccentricity (4.2.3) - - int i = 1; - while( i <= order ) - { - E = NKep( E, M, ecc ); - i *= 2; - } - - cout << "Exercise 4.2.2: Expansion of the Anomaly" << endl - << "Resulting expansion:" << endl << E << endl - << "Residual error:" << endl << (E-ecc*sin(E)-M) << endl; - - // sample the resulting polynomial over M0+-3 rad - ofstream f( "ex4_2_2.dat" ); - for( i = -300; i < 300; i++ ) - { - f << M0+i/100.0 << " " << E.evalScalar( i/100.0 ) - << " " << Kepler( M0+i/100.0, ecc0 ) << endl; - } - f.close( ); - // gnuplot command: plot 'ex4_2_2.dat' u ($1*180/pi):($2*180/pi) w l t 'DA', 'ex4_2_2.dat' u ($1*180/pi):($3*180/pi) w l t 'pointwise' -} - -int main( void ) -{ - DA::init( order, 2 ); // init with maximum computation order - - ex4_1_3( ); - ex4_1_4( ); - - ex4_2_1( 9.0 ); - ex4_2_2( 0.0, 0.5 ); -} \ No newline at end of file diff --git a/Tutorials/Tutorial2/5Vectors-Ex.cpp b/Tutorials/Tutorial2/5Vectors-Ex.cpp deleted file mode 100644 index 595f37b..0000000 --- a/Tutorials/Tutorial2/5Vectors-Ex.cpp +++ /dev/null @@ -1,111 +0,0 @@ -#include -#include -#include -using namespace std; using namespace DACE; - -const unsigned int order = 20; - -// Exercise 5.1.1: tangents and normals -template T f( AlgebraicVector x ) -{ - return -1/3*(sqr(x[0])+sqr(x[1])/2)+exp(x[0]/2+x[1]/4); -} - -void ex5_1_1( ) -{ - AlgebraicVector surf(3); - surf[0] = DA(1); - surf[1] = DA(2); - surf[2] = f( surf ); // trick: f() only uses the first 2 components of surf, which we already set - AlgebraicVector t1 = surf.deriv(1); - AlgebraicVector t2 = surf.deriv(2); - AlgebraicVector n = t1.cross(t2).normalize(); // normalized surface normal - - cout << "Exercise 5.1.1: tangents and normals" << endl - << t1 << t2 << n; -} - -// Exercise 5.1.2: (Uncontrolled) Equations of motion of the inverted pendulum -template AlgebraicVector ode_pendulum( AlgebraicVector x/*, T u*/ ) -{ - // constants - static const double l = 1.0; // length of pendulum (m) - static const double m = 0.1; // weight of balanced object (kg) - static const double M = 0.4; // weight of cart (kg) - static const double g = 9.81; // gravity acceleration constant on earth (kg*m/s^2) - - // variables - AlgebraicVector rhs(2); // right hand side of the ODE (2 dimensional) - T sint = sin(x[0]); // sine of theta - T cost = cos(x[0]); // cosine of theta - - // Equations of motion - rhs[0] = x[1]; - rhs[1] = (/*u+*/(M+m)*g*sint-m*l*sqr(x[1])*sint*cost)/((M+m)*l+m*l*sqr(cost)); - - return rhs; -} - -void ex5_1_2( ) -{ - AlgebraicVector x(2); - x[0] = 1.0; x[1] = 0.0; - cout << "Exercise 5.1.2: Equations of Motion" << endl << ode_pendulum( x ); -} - -// Exercise 5.2.1: Solar flux -void ex5_2_1( ) -{ - AlgebraicVector surf(3); - surf[0] = DA(1); - surf[1] = DA(2); - surf[2] = f( surf ); // trick: f() only uses the first 2 components of surf, which we already set - AlgebraicVector t1 = surf.deriv(1).normalize(); // normalizing these helps keep the coefficents small and prevents roundoff errors - AlgebraicVector t2 = surf.deriv(2).normalize(); - AlgebraicVector n = t1.cross(t2).normalize(); // normalized surface normal - AlgebraicVector sun(3); // sun direction - sun[0] = 0.0; sun[1] = 0.0; sun[2] = 1.0; - DA flux = n.dot(sun); // solar flux on the surface - - // Output results - cout << "Exercise 5.2.1: Solar flux" << endl << flux; - const int N = 30; - AlgebraicVector arg(2); - AlgebraicVector res(3); - ofstream file("ex5_2_1.dat"); - for( int i = -N; i <= N; i++ ) - { - arg[0] = (double)i/N; - for( int j = -N; j <= N; j++ ) - { - arg[1] = (double)j/N; - res = surf.eval(arg); - file << res[0] << " " << res[1] << " " << f( res ) << " " << flux.eval( arg ) << endl; - } - file << endl; - } - file.close( ); -} - -// Exercise 5.2.2: Area -void ex5_2_2( ) -{ - DA x = DA(1); - DA t = DA(2); - AlgebraicVector res(2); - res[0] = DA(1); - res[1] = ((1.0-x*x)*(t+1.0)+(x*x*x-x)*(1.0-t))/2.0; - - cout << "Exercise 5.2.2: Area" << endl << res; -} - -int main( void ) -{ - DA::init( order, 2 ); // init with maximum computation order - - ex5_1_1( ); - ex5_1_2( ); - - ex5_2_1( ); - ex5_2_2( ); -} diff --git a/Tutorials/Tutorial2/5Vectors-Ex.plot b/Tutorials/Tutorial2/5Vectors-Ex.plot deleted file mode 100644 index c004157..0000000 --- a/Tutorials/Tutorial2/5Vectors-Ex.plot +++ /dev/null @@ -1,11 +0,0 @@ -set term pdf -set output 'ex5_2_1.pdf' -#set term x11 - -set pm3d -set hidden3d -set view equal xyz -set title 'Solar flux' -splot 'ex5_2_1.dat' u 1:2:3:4 w l not - -#pause -1 diff --git a/Tutorials/Tutorial2/6Integrator-Ex.cpp b/Tutorials/Tutorial2/6Integrator-Ex.cpp deleted file mode 100644 index e89b2af..0000000 --- a/Tutorials/Tutorial2/6Integrator-Ex.cpp +++ /dev/null @@ -1,351 +0,0 @@ -#include -#include -#include -using namespace std; using namespace DACE; - -// Exercise 6.1.1: The Mid-point rule integrator -template T midpoint( T x0, double t0, double t1, T (*f)(T,double) ) -{ - const double hmax = 0.005; - int steps = ceil( (t1-t0)/hmax ); - double h = (t1-t0)/steps; - double t = t0; - - T xmid; - for( int i = 0; i < steps; i++ ) - { - xmid = x0 + 0.5*h*f( x0, t ); - x0 = x0 + h*f( xmid, t+0.5*h ); - t += h; - } - - return x0; -} - - -// Exercise 6.1.2: Model of the (uncontrolled) pendulum -// x = ( theta, theta_dot ) -// since the motion for x decouples we ignore it here -template AlgebraicVector pendulumRHS( AlgebraicVector x, double t ) -{ - // pendulum constants - static const double l = 0.61; // m - static const double m = 0.21; // kg - static const double M = 0.4926; // kg - static const double g = 9.81; // kg*m/s^2 - - AlgebraicVector res(2); - - res[0] = x[1]; - res[1] = ((M+m)*g*sin(x[0])-m*l*sqr(x[1])*sin(x[0])*cos(x[0]))/((M+m)*l+m*l*sqr(cos(x[0]))); - - return res; -} - -void ex6_1_2( ) -{ - AlgebraicVector xDA(2); - AlgebraicVector xdb(2); - double t; - const double dt = 0.05; // take a snap shot every 0.1s - - ofstream file( "ex6_1_2.dat" ); - file.precision( 16 ); - - xdb[0] = 0.0; xdb[1] = 0.2; // initial condition (double) - t = 0; - for( int i = 0; i < 100; i++ ) - { - xdb = midpoint( xdb, t, t+dt, pendulumRHS ); // propagate forward for dt seconds - file << t << " " << xdb[0] << " " << xdb[1] << endl; - t += dt; - } - file << endl << endl; - - xDA[0] = 0.0; xDA[1] = 0.2+0.04*DA(1); // initial condition (DA) - t = 0; - for( int i = 0; i < 100; i++ ) - { - xDA = midpoint( xDA, t, t+dt, pendulumRHS ); // propagate forward for dt seconds - file << t << " " << xDA[0].evalScalar(-1.0) << " " << xDA[0].evalScalar(+1.0); - file << " " << xDA[1].evalScalar(-1.0) << " " << xDA[1].evalScalar(+1.0) << endl; - t += dt; - } - file << endl; - file.close( ); - - cout << "Exercise 6.1.2: Model of the (uncontrolled) pendulum" << endl << xDA << endl; -} - - -// Exercise 6.1.3: Set propagation -// the right hand side -template AlgebraicVector f( AlgebraicVector x, double t ) -{ - const double alpha = 0.1; - AlgebraicVector res(2); - - res[0] = -x[1]; - res[1] = x[0]; - return (1.0 + alpha*x.vnorm())*res; -} - -// convenience routine to evaluate and plot -void plot( AlgebraicVector x, double t, int N, ofstream &file ) -{ - AlgebraicVector arg(2), res; - for( int i = -N; i <= N; i++ ) - { - arg[0] = (double)i/N; - for( int j = -N; j <= N; j++ ) - { - arg[1] = (double)j/N; - res = x.eval(arg); - file << t << " " << res[0] << " " << res[1] << endl; - } - file << endl; - } - file << endl; -} - -void ex6_1_3( ) -{ - const double pi = 3.141592653589793; - AlgebraicVector x(2); - double t; - const double dt = 2.0*pi/6.0; - - ofstream file( "ex6_1_3.dat" ); - file.precision( 16 ); - - x[0] = 2.0 + DA(1); x[1] = DA(2); // initial condition box - t = 0; - plot( x, t, 7, file ); - - for( int i = 0; i < 6; i++ ) - { - x = midpoint( x, t, t+dt, f ); // propagate forward for dt seconds - t += dt; - plot( x, t, 7, file ); - } - file.close( ); - - cout << "Exercise 6.1.3: Set propagation" << endl << x << endl; -} - - -// Exercise 6.1.4: State Transition Matrix -void ex6_1_4( ) -{ - const double pi = 3.141592653589793; - AlgebraicVector x(2); - - x[0] = 1.0 + DA(1); x[1] = 1.0 + DA(2); // initial condition around (1,1) - x = midpoint( x, 0, 2*pi, f ); - - cout << "Exercise 6.1.4: State Transition Matrix" << endl; - cout.precision( 7 ); - cout << x[0].deriv(1).cons() << " " << x[0].deriv(2).cons() << endl; - cout << x[1].deriv(1).cons() << " " << x[1].deriv(2).cons() << endl; - cout << endl; -} - - -// Exercise 6.1.5: Parameter dependence -// the right hand side (note: now it can only be evaluated with DA because alpha is a DA!) -AlgebraicVector fParam( AlgebraicVector x, double t ) -{ - const DA alpha = 0.05 + 0.05*DA(1); // parameter, now it's a DA - AlgebraicVector res(2); - - res[0] = -x[1]; - res[1] = x[0]; - return (1.0 + alpha*x.vnorm())*res; -} - -void ex6_1_5( ) -{ - const double pi = 3.141592653589793; - AlgebraicVector x(2); - - x[0] = 1.0; x[1] = 1.0; // initial condition (1,1) - x = midpoint( x, 0, 2*pi, fParam ); - - ofstream file( "ex6_1_5.dat" ); - file.precision( 16 ); - file << "1 1" << endl << endl << endl; - for( int i = 0; i <= 20; i++ ) - file << x[0].evalScalar(-1.0+i/10.0) << " " << x[1].evalScalar(-1.0+i/10.0) << endl; - file.close( ); - - cout << "Exercise 6.1.5: Parameter dependence" << endl << x; - cout << endl; -} - - -// Exercise 6.2.1: 3/8 rule RK4 integrator -template T rk4( T x0, double t0, double t1, T (*f)(T,double) ) -{ - const double hmax = 0.01; - int steps = ceil( (t1-t0)/hmax ); - double h = (t1-t0)/steps; - double t = t0; - - T k1, k2, k3, k4; - for( int i = 0; i < steps; i++ ) - { - k1 = f( x0, t ); - k2 = f( x0 + h*k1/3.0, t + h/3.0 ); - k3 = f( x0 + h*(-k1/3.0 + k2), t + 2.0*h/3.0 ); - k4 = f( x0 + h*(k1 - k2 + k3), t + h ); - x0 = x0 + h*(k1 + 3*k2 + 3*k3 +k4)/8.0; - t += h; - } - - return x0; -} - - -// Exercise 6.2.2: Artsy Set Propagation -void ex6_2_2( ) -{ - const double pi = 3.141592653589793; - AlgebraicVector x(2); - double t; - const double dt = 2.0*pi/6.0; - - ofstream file( "ex6_2_2.dat" ); - file.precision( 16 ); - - // initial condition (c.f. 5Vectors-Ex.cpp) - x[0] = DA(1); - x[1] = ((1.0-DA(1)*DA(1))*(DA(2)+1.0)+(DA(1)*DA(1)*DA(1)-DA(1))*(1.0-DA(2)))/2.0; - t = 0; - plot( x, t, 7, file ); - - for( int i = 0; i < 6; i++ ) - { - x = midpoint( x, t, t+dt, f ); // propagate forward for dt seconds - t += dt; - plot( x, t, 7, file ); - } - file.close( ); - - cout << "Exercise 6.2.2: Artsy Set propagation" << endl << x << endl; -} - - -// Exercise 6.2.3: CR3BP -template AlgebraicVector CR3BP( AlgebraicVector x, double t ) -{ - const double MU = 0.30404234e-5; - AlgebraicVector res(6); - T d1, d2; - - d1 = sqrt( sqr(x[0]+MU) + sqr(x[1]) + sqr(x[2]) ); - d1 = 1.0/(d1*d1*d1); // first distance - d2 = sqrt( sqr(x[0]+MU-1.0) + sqr(x[1]) + sqr(x[2]) ); - d2 = 1.0/(d2*d2*d2); // second distance - - res[0] = x[3]; - res[1] = x[4]; - res[2] = x[5]; - res[3] = x[0] + 2.0*x[4] - d1*(1-MU)*(x[0]+MU) - d2*MU*(x[0]+MU-1.0); - res[4] = x[1] - 2.0*x[3] - d1*(1-MU)*x[1] - d2*MU*x[1]; - res[5] = - d1*(1-MU)*x[2] - d2*MU*x[2]; - - return res; -} - -void ex6_2_3( ) -{ - double T = 3.05923; - AlgebraicVector x0(6); - x0[0] = 0.9888426847; x0[1] = 0; x0[2] = 0.0011210277; x0[3] = 0; x0[4] = 0.0090335498; x0[5] = 0; - AlgebraicVector x = x0 + AlgebraicVector::identity( ); - - DA::pushTO( 1 ); // only first order computation needed - - x = rk4( x, 0, T, CR3BP ); - - cout << "Exercise 6.2.3: CR3BP STM" << endl; - cout.precision( 6 ); - cout << cons(x); - for( int i = 0; i < 6; i++ ) - { - for( int j = 1; j <= 6; j++ ) - { - cout << cons(x[i].deriv(j)) << " "; - } - cout << endl; - } - cout << endl; - - DA::popTO( ); -} - - -// Exercise 6.2.4: Set propagation revisited -void ex6_2_4( ) -{ - const double pi = 3.141592653589793; - AlgebraicVector x(2); - double t; - const double dt = 2.0*pi/6.0; - - ofstream file( "ex6_2_4.dat" ); - file.precision( 16 ); - - // initial condition box, in polar coordinates - x[0] = cos(0.3*DA(2))*(2.0 + DA(1)); x[1] = sin(0.3*DA(2))*(2.0 + DA(1)); - t = 0; - plot( x, t, 40, file ); - - for( int i = 0; i < 6; i++ ) - { - x = midpoint( x, t, t+dt, f ); // propagate forward for dt seconds - t += dt; - plot( x, t, 40, file ); - } - file.close( ); - - cout << "Exercise 6.2.4: Set propagation revisited" << endl << x << endl; -} - - -// Exercise 6.2.5: The State Transition Matrix reloaded -void ex6_2_5( ) -{ - const double pi = 3.141592653589793; - AlgebraicVector x(2); - - x[0] = 1.0 + DA(2); x[1] = 1.0 + DA(3); // initial condition (1,1) plus DA identity (but in DA(2) and DA(3) as DA(1) is already used for alpha!) - x = midpoint( x, 0, 2*pi, fParam ); - - AlgebraicVector arg(3); - arg[0] = DA(1); // we want to evaluate the derivatives at (alpha,0,0), so keep DA(1) and replace DA(2) and DA(3) by zero - arg[1] = 0; - arg[2] = 0; - - cout << "Exercise 6.2.5: The State Transition Matrix reloaded" << endl; - cout << x[0].deriv(2).eval(arg) << x[0].deriv(3).eval(arg) << endl; - cout << x[1].deriv(2).eval(arg) << x[1].deriv(3).eval(arg) << endl; - cout << endl; -} - - -int main( void ) -{ - DA::init( 15, 6 ); // init with maximum computation order - - ex6_1_2( ); - ex6_1_3( ); - ex6_1_4( ); - ex6_1_5( ); - - ex6_2_2( ); - ex6_2_3( ); - ex6_2_4( ); - ex6_2_5( ); -} - diff --git a/Tutorials/Tutorial2/6Integrator-Ex.plot b/Tutorials/Tutorial2/6Integrator-Ex.plot deleted file mode 100644 index 7844343..0000000 --- a/Tutorials/Tutorial2/6Integrator-Ex.plot +++ /dev/null @@ -1,61 +0,0 @@ -set term pdf -set output 'ex6_1_2.pdf' - -set xlabel 'time [s]' -set ylabel 'angle [rad]' -plot 'ex6_1_2.dat' using 1:2 i 0 with lines notitle -set ylabel 'angular velocity [rad/s]' -plot 'ex6_1_2.dat' using 1:3 i 0 with lines notitle -set ylabel 'angle [rad]' -plot 'ex6_1_2.dat' using 1:2 i 1 with lines title 'left boundary', 'ex6_1_2.dat' using 1:3 i 1 with lines title 'right boundary' -set ylabel 'angular velocity [rad/s]' -plot 'ex6_1_2.dat' using 1:4 i 1 with lines title 'left boundary', 'ex6_1_2.dat' using 1:5 i 1 with lines title 'right boundary' - - -set term pdf -set output 'ex6_1_3.pdf' - -unset xlabel -unset ylabel -set key bmargin -set view 0, 0 -set view equal xyz -splot for [i=0:6] 'ex6_1_3.dat' using 2:3:(0) index i with lines title 'snapshot '.i - - -set term pdf -set output 'ex6_1_5.pdf' - -unset xlabel -unset ylabel -set size square -set key rmargin -plot [-2:2] [-2:2] \ - 'ex6_1_5.dat' index 0 title 'Starting point', \ - 'ex6_1_5.dat' index 1 with line title 'Final states' - - -set term pdf -set output 'ex6_2_2.pdf' - -unset xlabel -unset ylabel -set key bmargin -set view 0, 0 -set view equal xyz -splot [-2:2] [-2:2] for [i=0:6] 'ex6_2_2.dat' using 2:3:(0) index i with lines title 'snapshot '.i - -do for [i=0:6] { - splot [-2:2] [-2:2] 'ex6_2_2.dat' using 2:3:(0) index i with lines title 'snapshot '.i -} - - -set term pdf -set output 'ex6_2_4.pdf' - -unset xlabel -unset ylabel -set key bmargin -set view 0, 0 -set view equal xyz -splot for [i=0:6] 'ex6_2_4.dat' using 2:3:(0) index i with lines title 'snapshot '.i diff --git a/Tutorials/Tutorial2/7PID-Ex.cpp b/Tutorials/Tutorial2/7PID-Ex.cpp deleted file mode 100644 index 3f9db4b..0000000 --- a/Tutorials/Tutorial2/7PID-Ex.cpp +++ /dev/null @@ -1,265 +0,0 @@ -#include -#include -#include -#include -using namespace DACE; using namespace std; - -// Exercise 6.2.1: 3/8 rule RK4 integrator -template T rk4( T x0, double t0, double t1, T (*f)(T,double) ) -{ - const double hmax = 0.01; - int steps = ceil( (t1-t0)/hmax ); - double h = (t1-t0)/steps; - double t = t0; - - T k1, k2, k3, k4; - for( int i = 0; i < steps; i++ ) - { - k1 = f( x0, t ); - k2 = f( x0 + h*k1/3.0, t + h/3.0 ); - k3 = f( x0 + h*(-k1/3.0 + k2), t + 2.0*h/3.0 ); - k4 = f( x0 + h*(k1 - k2 + k3), t + h ); - x0 = x0 + h*(k1 + 3*k2 + 3*k3 +k4)/8.0; - t += h; - } - - return x0; -} - - -// Exercise 7.1.1: Model of the controlled pendulum -// x = ( theta, theta_dot, u ) -template AlgebraicVector pendulumRHS( AlgebraicVector x, double t ) -{ - // pendulum constants - static const double l = 0.61; // m - static const double m = 0.21; // kg - static const double M = 0.4926; // kg - static const double g = 9.81; // kg*m/s^2 - - AlgebraicVector res(3); - T sint, cost; - sint = sin(x[0]); - cost = cos(x[0]); - - res[0] = x[1]; - res[1] = (x[2]+(M+m)*g*sint-m*l*sqr(x[1])*sint*cost)/((M+m)*l+m*l*sqr(cost)); - res[2] = 0; // u is assumed constant unless changed externally by the controller - - return res; -} - - -// Exercise 7.1.2: Tuning the PID simulator (double) -void ex7_1_2( ) -{ - const double Kp = 8, Ti = 3, Td = 0.3; // PID parameters - const double setPt = 0.0, dt = 0.05; // Set point, controller time step (50 ms) - const double umax = 10; // maximum control (Exercise 7.2.1) - - ofstream file( "ex7_1_2.dat" ); - - double intErr = 0, lastErr = 0; // PID controller variables - - double t = 0; - AlgebraicVector x(3); - x[0] = 0.1; x[1] = 0; x[2] = 0; // Initial condition (u(0)=x[2]=0) - - // propagate the model for 100 sec - while( (t<100.0) && (abs(x[0])<1.5) ) - { - // compute the PID control at this controller time step - double err = setPt-x[0]; - double derr = (err-lastErr)/dt; - intErr += lastErr*dt; - lastErr = err; - x[2] = Kp*(err + Td*derr + intErr/Ti); - // prevent control saturation (Exercise 7.2.1) -// x[2] = tanh(x[2]/2.0/umax)*umax*2.0; -// x[2] = max(min(x[2],umax),-umax); - - // output and propagate one time step - file << t << " " << x[0] << " " << x[2] << endl; - x = rk4( x, t, t+dt, pendulumRHS ); - t += dt; - } - - cout << "Final angle theta:" << x << endl; - if(abs(x[0])>1.5) - cout << "WHOOPSY: Fell over after " << t << " seconds." << endl; - - file.close(); -} - - -// Exercise 7.1.3: PID simulator (DA) -void ex7_1_3( ) -{ - const double Kp = 8, Ti = 3, Td = 0.3; // PID parameters - const double setPt = 0.0, dt = 0.05; // Set point, controller time step (50 ms) - const double umax = 10; // maximum control (Exercise 7.2.1) - - ofstream file( "ex7_1_3.dat" ); - - DA intErr = 0, lastErr = 0; // PID controller variables - - double t = 0; - AlgebraicVector x(3); - x[0] = 0.1+0.1*DA(1); x[1] = 0; x[2] = 0; // Initial condition - - // propagate the model state for 100 sec - while( (t<40.0) && (abs(cons(x[0]))<1.5) ) - { - // compute the PID control - DA err = setPt-x[0]; - DA derr = (err-lastErr)/dt; - intErr += lastErr*dt; - lastErr = err; - x[2] = Kp*(err + Td*derr + intErr/Ti); - // prevent control saturation (Exercise 7.2.1) -// x[2] = tanh(x[2]/umax)*umax; - - // output and propagate one time step (Exercise 7.1.4) - Interval bx = x[0].bound( ); - Interval bu = x[2].bound( ); - file << t - << " " << cons(x[0]) << " " << bx.m_ub << " " << bx.m_lb - << " " << x[0].evalScalar(-1.0) << " " << x[0].evalScalar(1.0) - << " " << cons(x[2]) << " " << bu.m_ub << " " << bu.m_lb - << " " << x[2].evalScalar(-1.0) << " " << x[2].evalScalar(1.0) << endl; - x = rk4( x, t, t+dt, pendulumRHS ); - t += dt; - } - - cout << "Final angle theta:" << x << endl; - if(abs(cons(x[0]))>1.5) - cout << "WHOOPSY: Fell over after " << t << " seconds." << endl; - - file.close(); -} - - -// Exercise 7.1.5: Bounding -void ex7_1_5( ) -{ - DA x = DA(1); DA y = DA(2); - DA func = sin(x/2)/(2+cos(y/2+x*x)); - - Interval a, b, c; - - // bound by rasterizing - AlgebraicVector arg(2); - a.m_lb = 9999999; a.m_ub = -a.m_lb; - c = a; - for( int i = -10; i < 10; i++ ) - { - arg[0] = i/10.0; - for( int j = -10; j < 10; j++ ) - { - arg[1] = j/10.0; - // polynomial expansion - double r = func.eval( arg ); - a.m_lb = min( a.m_lb, r ); - a.m_ub = max( a.m_ub, r ); - // actual function - r = sin(arg[0]/2)/(2+cos(arg[1]/2+arg[0]*arg[0])); - c.m_lb = min( c.m_lb, r ); - c.m_ub = max( c.m_ub, r ); - } - } - - // DA bounding - b = func.bound( ); - - cout.precision( 16 ); - cout << "func:" << endl << func; - cout << "Bounds:" << endl - << "DA bound: [" << b.m_lb << "," << b.m_ub << "]" << endl - << "DA raster: [" << a.m_lb << "," << a.m_ub << "]" << endl - << "double raster: [" << c.m_lb << "," << c.m_ub << "]" << endl << endl; -} - - -// Exercise 7.2.2: PID simulator with uncertain mass (DA) - -// Model of controlled pendulum with uncertain mass -// x = ( theta, theta_dot, u ) -AlgebraicVector pendulumRHSmass( AlgebraicVector x, double t ) -{ - // pendulum constants - static const double l = 0.61; // m - static const DA m = 0.21*(1+0.1*DA(1)); // kg - static const double M = 0.4926; // kg - static const double g = 9.81; // kg*m/s^2 - - AlgebraicVector res(3); - DA sint, cost; - sint = sin(x[0]); - cost = cos(x[0]); - - res[0] = x[1]; - res[1] = (x[2]+(M+m)*g*sint-m*l*sqr(x[1])*sint*cost)/((M+m)*l+m*l*sqr(cost)); - res[2] = 0; // u is assumed constant unless changed externally by the controller - - return res; -} - -void ex7_2_2( ) -{ - const double Kp = 8, Ti = 3, Td = 0.3; // PID parameters - const double setPt = 0.0, dt = 0.05; // Set point, controller time step (50 ms) - const double umax = 10; // maximum control (Exercise 7.2.1) - - ofstream file( "ex7_2_2.dat" ); - - DA intErr = 0, lastErr = 0; // PID controller variables - - double t = 0; - AlgebraicVector x(3); - x[0] = 0.1; x[1] = 0; x[2] = 0; // Initial condition - - // propagate the model state for 100 sec - while( (t<40.0) && (abs(cons(x[0]))<1.5) ) - { - // compute the PID control - DA err = setPt-x[0]; - DA derr = (err-lastErr)/dt; - intErr += lastErr*dt; - lastErr = err; - x[2] = Kp*(err + Td*derr + intErr/Ti); - // prevent control saturation (Exercise 7.2.1) -// x[2] = tanh(x[2]/umax)*umax; - - // output and propagate one time step (Exercise 7.1.4) - Interval bx = x[0].bound( ); - Interval bu = x[2].bound( ); - file << t - << " " << cons(x[0]) << " " << bx.m_ub << " " << bx.m_lb - << " " << x[0].evalScalar(-1.0) << " " << x[0].evalScalar(1.0) - << " " << cons(x[2]) << " " << bu.m_ub << " " << bu.m_lb - << " " << x[2].evalScalar(-1.0) << " " << x[2].evalScalar(1.0) << endl; - x = rk4( x, t, t+dt, pendulumRHSmass ); - t += dt; - } - - cout << "Final angle theta:" << x << endl; - if(abs(cons(x[0]))>1.5) - cout << "WHOOPSY: Fell over after " << t << " seconds." << endl; - - file.close(); -} - - - -int main( int argc, char** argv ) -{ - DA::init( 10, 2 ); - - ex7_1_2( ); - ex7_1_3( ); - ex7_1_5( ); - - ex7_2_2( ); - - return 0; -} diff --git a/Tutorials/Tutorial2/7PID-Ex.plot b/Tutorials/Tutorial2/7PID-Ex.plot deleted file mode 100644 index 538842a..0000000 --- a/Tutorials/Tutorial2/7PID-Ex.plot +++ /dev/null @@ -1,46 +0,0 @@ -deg(x)=x*180/pi - - -set term pdf -set output 'ex7_1_2.pdf' - -set xlabel 'time [s]' -set ylabel 'angle [rad]' -plot 'ex7_1_2.dat' u 1:(deg($2)) w l t 'angle' -set ylabel 'force [N]' -plot 'ex7_1_2.dat' u 1:3 w l t 'control' - - -set term pdf -set output 'ex7_1_3.pdf' - -set ylabel 'angle [rad]' -plot \ - 'ex7_1_3.dat' u 1:(deg($3)):(deg($4)) w filledcurves lc 'gray' not,\ - 'ex7_1_3.dat' u 1:(deg($5)) w l t 'angle (left boundary)',\ - 'ex7_1_3.dat' u 1:(deg($6)) w l t 'angle (right boundary)',\ - 'ex7_1_3.dat' u 1:(deg($2)) w l t 'angle (reference)' -set ylabel 'force [N]' -plot \ - 'ex7_1_3.dat' u 1:8:9 w filledcurves lc 'gray' not,\ - 'ex7_1_3.dat' u 1:10 w l t 'control (left boundary)',\ - 'ex7_1_3.dat' u 1:11 w l t 'control (right boundary)',\ - 'ex7_1_3.dat' u 1:7 w l t 'control (reference)' - - -set term pdf -set output 'ex7_2_2.pdf' - -set ylabel 'angle [rad]' -plot \ - 'ex7_2_2.dat' u 1:(deg($3)):(deg($4)) w filledcurves lc 'gray' not,\ - 'ex7_2_2.dat' u 1:(deg($5)) w l t 'angle (left boundary)',\ - 'ex7_2_2.dat' u 1:(deg($6)) w l t 'angle (right boundary)',\ - 'ex7_2_2.dat' u 1:(deg($2)) w l t 'angle (reference)' -set ylabel 'force [N]' -plot \ - 'ex7_2_2.dat' u 1:8:9 w filledcurves lc 'gray' not,\ - 'ex7_2_2.dat' u 1:10 w l t 'control (left boundary)',\ - 'ex7_2_2.dat' u 1:11 w l t 'control (right boundary)',\ - 'ex7_2_2.dat' u 1:7 w l t 'control (reference)' - diff --git a/Tutorials/Tutorial2/CMakeLists.txt b/Tutorials/Tutorial2/CMakeLists.txt deleted file mode 100644 index bab3247..0000000 --- a/Tutorials/Tutorial2/CMakeLists.txt +++ /dev/null @@ -1,49 +0,0 @@ -# DACE Examples -# -# This file shows how to use the DACE in a cmake project -# Simply import the DACE cmake package. It automatically provides the correct -# locations for headers and libraries. -# If the DACE is installed in a non-default location, tell cmake where to look: -# set(dace_DIR "/my/custom/root/usr/local/lib/cmake/dace") -# It exposes two targets: -# dace::dace (the dynamic DACE library) -# dace::dace_s (the static DACE library) -# Make sure to link your executable with one of those and include the correct -# header ( or ) in the code. -# -# On Windows only, the dace.dll dynamics library is typically expected to be -# located in the same directory as the executable using it. You can use a -# command such as the file(COPY ...) below to automatically copy the dace.dll -# from the DACE package next to your executables. - -cmake_minimum_required (VERSION 2.8.4) - -project(Examples2 CXX) - -find_package(dace 2.0.0 REQUIRED) - -if(WIN32) - get_target_property(DACEDLL dace::dace LOCATION) - file(COPY "${DACEDLL}" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}") -endif(WIN32) - -add_executable(1Basics-Ex 1Basics-Ex.cpp) -target_link_libraries(1Basics-Ex PUBLIC dace::dace) - -add_executable(2IntDiff-Ex 2IntDiff-Ex.cpp) -target_link_libraries(2IntDiff-Ex PUBLIC dace::dace) - -add_executable(3Evaluation-Ex 3Evaluation-Ex.cpp) -target_link_libraries(3Evaluation-Ex PUBLIC dace::dace) - -add_executable(4Solver-Ex 4Solver-Ex.cpp) -target_link_libraries(4Solver-Ex PUBLIC dace::dace) - -add_executable(5Vectors-Ex 5Vectors-Ex.cpp) -target_link_libraries(5Vectors-Ex PUBLIC dace::dace) - -add_executable(6Integrator-Ex 6Integrator-Ex.cpp) -target_link_libraries(6Integrator-Ex PUBLIC dace::dace) - -add_executable(7PID-Ex 7PID-Ex.cpp) -target_link_libraries(7PID-Ex PUBLIC dace::dace) diff --git a/Tutorials/Tutorial2/Tutorial2.pdf b/Tutorials/Tutorial2/Tutorial2.pdf deleted file mode 100644 index 856c089..0000000 Binary files a/Tutorials/Tutorial2/Tutorial2.pdf and /dev/null differ diff --git a/core/include/dace/daceerror.h b/core/include/dace/daceerror.h index e227d44..a28d140 100644 --- a/core/include/dace/daceerror.h +++ b/core/include/dace/daceerror.h @@ -48,7 +48,6 @@ DACE_API const errstrings DACEerr[] = { {1006, "Incorrect DA coding arrays"}, {1007, "Requested length too long"}, {1008, "Error in monomial evaluation tree construction"}, - { 8, ""}, { 9, ""}, { 10, ""}, { 911, "Order and/or variable too large"},