source: trunk/test_more/div_eq.cpp @ 2506

Last change on this file since 2506 was 2506, checked in by bradbell, 8 years ago

Change Licenses: CPL-1.0 -> EPL-1.0, GPL-2.0->GPL-3.0

  • Property svn:keywords set to Id
File size: 3.9 KB
Line 
1/* $Id: div_eq.cpp 2506 2012-10-24 19:36:49Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse 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/*
14Old DivEq example now used just for valiadation testing
15*/
16# include <cppad/cppad.hpp>
17
18namespace { // BEGIN empty namespace
19
20bool DivEqTestOne(void)
21{       bool ok = true;
22
23        using namespace CppAD;
24
25        // independent variable vector, indices, values, and declaration
26        CPPAD_TESTVECTOR(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_TESTVECTOR(AD<double>) Z(2);
35        size_t x = 0;
36        size_t y = 1;
37
38        // some constants
39        AD<double> zero = 0.;
40        AD<double> one  = 1.;
41
42        // dependent variable values
43        Z[x] = U[s];   
44        Z[y] = U[t];
45        Z[x] /= U[t]; // AD<double> *= AD<double>
46        Z[y] /= 5.;   // AD<double> *= double
47        zero /= Z[y]; // divide into a parameter equal to zero
48        Z[y] /= one;  // divide by a parameter equal to one
49        Z[y] /= 1.;   // divide by a double equal to one
50       
51        // check that zero is still a parameter
52        // (must do this before creating f because it erases the tape)
53        ok &= Parameter(zero);
54
55        // create f : U -> Z and vectors for derivative calcualtions
56        ADFun<double> f(U, Z);
57        CPPAD_TESTVECTOR(double) v( f.Domain() );
58        CPPAD_TESTVECTOR(double) w( f.Range() );
59
60        // check that none of the components of f are parameters
61        size_t i;
62        for(i = 0; i < f.Range(); i++)
63                ok &= ! f.Parameter(i);
64
65        // check functin values
66        ok &= NearEqual(Z[x] , 3. / 2. , 1e-10, 1e-10);
67        ok &= NearEqual(Z[y] , 2. / 5. , 1e-10, 1e-10);
68
69        // forward computation of partials w.r.t. t
70        v[s] = 0.;
71        v[t] = 1.;
72        w    = f.Forward(1, v);
73        ok &= NearEqual(w[x] , -1.*U[s]/(U[t]*U[t]) , 1e-10, 1e-10); // dx/dt
74        ok &= NearEqual(w[y] , 1. / 5.              , 1e-10, 1e-10); // dy/dt
75
76        // reverse computation of second partials of x
77        CPPAD_TESTVECTOR(double) r( f.Domain() * 2 );
78        w[x] = 1.;
79        w[y] = 0.;
80        r    = f.Reverse(2, w); 
81        ok &= NearEqual(r[2 * s + 1]                 // d^2 x / (dt ds)
82                 , - 1. / (U[t] * U[t])     , 1e-10 , 1e-10 );
83        ok &= NearEqual(r[2 * t + 1]                 // d^2 x / (dt dt)
84                , 2. * U[s] / (U[t] * U[t] * U[t])  , 1e-10 , 1e-10 );
85
86        return ok;
87}
88
89bool DivEqTestTwo(void)
90{       bool ok = true;
91
92        using namespace CppAD;
93
94        // independent variable vector
95        double u0 = .5;
96        CPPAD_TESTVECTOR(AD<double>) U(1);
97        U[0]      = u0; 
98        Independent(U);
99
100        // dependent variable vector
101        CPPAD_TESTVECTOR(AD<double>) Z(1);
102        Z[0] = U[0] * U[0]; // initial value
103        Z[0] /= 2;          // AD<double> /= int
104        Z[0] /= 4.;         // AD<double> /= double
105        Z[0] /= U[0];       // AD<double> /= AD<double>
106
107        // create f: U -> Z and vectors used for derivative calculations
108        ADFun<double> f(U, Z); 
109        CPPAD_TESTVECTOR(double) v(1);
110        CPPAD_TESTVECTOR(double) w(1);
111
112        // check value
113        ok &= NearEqual(Z[0] , u0*u0/(2*4*u0),  1e-10 , 1e-10);
114
115        // forward computation of partials w.r.t. u
116        size_t j;
117        size_t p     = 5;
118        double jfac  = 1.;
119        double value = 1./8.;
120        v[0]         = 1.;
121        for(j = 1; j < p; j++)
122        {       jfac *= j;
123                w     = f.Forward(j, v);       
124                ok &= NearEqual(jfac*w[0], value, 1e-10 , 1e-10); // d^jz/du^j
125                v[0]  = 0.;
126                value = 0.;
127        }
128
129        // reverse computation of partials of Taylor coefficients
130        CPPAD_TESTVECTOR(double) r(p); 
131        w[0]  = 1.;
132        r     = f.Reverse(p, w);
133        jfac  = 1.;
134        value = 1./8.;
135        for(j = 0; j < p; j++)
136        {       ok &= NearEqual(jfac*r[j], value, 1e-10 , 1e-10); // d^jz/du^j
137                jfac *= (j + 1);
138                value = 0.;
139        }
140
141        return ok;
142}
143
144} // END empty namespace
145
146bool DivEq(void)
147{       bool ok = true;
148        ok &= DivEqTestOne();
149        ok &= DivEqTestTwo(); 
150        return ok;
151}
Note: See TracBrowser for help on using the repository browser.