source: trunk/test_more/sub_eq.cpp @ 2354

Last change on this file since 2354 was 1370, checked in by bradbell, 11 years ago

trunk: Fix svn_add_id.sh and use it set Id property for some missed files.

  • Property svn:keywords set to Id
File size: 3.2 KB
Line 
1/* $Id: sub_eq.cpp 1370 2009-05-31 05:31:50Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Common Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13/*
14Two old SubEq examples now used just for valiadation testing
15*/
16# include <cppad/cppad.hpp>
17
18namespace { // BEGIN empty namespace
19
20bool SubEqTestOne(void)
21{       bool ok = true;
22
23        using namespace CppAD;
24
25        // independent variable vector, indices, values, and declaration
26        CPPAD_TEST_VECTOR< AD<double> > U(2);
27        size_t s = 0;
28        size_t t = 1;
29        U[s] = 3.;
30        U[t] = 2.;
31        Independent(U);
32
33        // dependent variable vector and indices
34        CPPAD_TEST_VECTOR< AD<double> > Z(2);
35        size_t x = 0;
36        size_t y = 1;
37
38        // dependent variable values
39        Z[x]  = U[s]; 
40        Z[y]  = U[t];
41        Z[x] -= U[t];  // AD<double> -= AD<double>
42        Z[y] -= 5.;    // AD<double> -= double
43
44        // create f: U -> Z and vectors used for derivative calculations
45        ADFun<double> f(U, Z);
46        CPPAD_TEST_VECTOR<double> v( f.Domain() );
47        CPPAD_TEST_VECTOR<double> w( f.Range() );
48
49        // check function values
50        ok &= ( Z[x] == 3. - 2. );
51        ok &= ( Z[y] == 2. - 5. );
52
53        // forward computation of partials w.r.t. t
54        v[s] = 0.;
55        v[t] = 1.;
56        w = f.Forward(1, v);
57        ok &= ( w[x] == -1. );  // dx/dt
58        ok &= ( w[y] == 1. );   // dy/dt
59
60        // reverse computation of second partials of x
61        CPPAD_TEST_VECTOR<double> r( f.Domain() * 2 );
62        w[x] = 1.;
63        w[y] = 0.;
64        r = f.Reverse(2, w);
65        ok &= ( r[2 * s + 1] == 0. );  // d^2 x / (ds ds)
66        ok &= ( r[2 * t + 1] == 0. );  // d^2 x / (ds dt)
67
68        return ok;
69}
70
71bool SubEqTestTwo(void)
72{       bool ok = true;
73
74        using namespace CppAD;
75
76        // independent variable vector
77        double u0 = .5;
78        CPPAD_TEST_VECTOR< AD<double> > U(1);
79        U[0]      = u0; 
80        Independent(U);
81
82        // dependent variable vector
83        CPPAD_TEST_VECTOR< AD<double> > Z(1);
84        Z[0] = U[0];       // initial value
85        Z[0] -= 2;         // AD<double> -= int
86        Z[0] -= 4.;        // AD<double> -= double
87        Z[0] -= 2 * U[0];  // AD<double> -= AD<double>
88
89        // create f: U -> Z and vectors used for derivative calculations
90        ADFun<double> f(U, Z); 
91        CPPAD_TEST_VECTOR<double> v(1);
92        CPPAD_TEST_VECTOR<double> w(1);
93
94        // check value
95        ok &= NearEqual(Z[0] , u0-2-4-2*u0,  1e-10 , 1e-10);
96
97        // forward computation of partials w.r.t. u
98        size_t j;
99        size_t p     = 5;
100        double jfac  = 1.;
101        double value = -1.;
102        v[0]         = 1.;
103        for(j = 1; j < p; j++)
104        {       jfac *= j;
105                w     = f.Forward(j, v);       
106                ok &= NearEqual(jfac*w[0], value, 1e-10 , 1e-10); // d^jz/du^j
107                v[0]  = 0.;
108                value = 0.;
109        }
110
111        // reverse computation of partials of Taylor coefficients
112        CPPAD_TEST_VECTOR<double> r(p); 
113        w[0]  = 1.;
114        r     = f.Reverse(p, w);
115        jfac  = 1.;
116        value = -1.;
117        for(j = 0; j < p; j++)
118        {       ok &= NearEqual(jfac*r[j], value, 1e-10 , 1e-10); // d^jz/du^j
119                jfac *= (j + 1);
120                value = 0.;
121        }
122
123        return ok;
124}
125
126} // END empty namespace
127
128bool SubEq(void)
129{       bool ok = true;
130        ok &= SubEqTestOne();
131        ok &= SubEqTestTwo(); 
132        return ok;
133}
Note: See TracBrowser for help on using the repository browser.