11
11
* -----------------------------------------------------------------------------
12
12
* This example solves the Lotka-Volterra ODE with four parameters,
13
13
*
14
- * u = [dx/dt] = [ p_0*x - p_1*x*y ]
14
+ * u = [dx/dt] = [ p_0*x - p_1*x*y ]
15
15
* [dy/dt] [ -p_2*y + p_3*x*y ].
16
16
*
17
- * The initial condition is u(t_0) = 1.0 and we use the parameters
18
- * p = [1.5, 1.0, 3.0, 1.0]. The integration interval is t \in [0, 10.].
19
- * The implicit BDF method from CVODES is used to solve the forward problem.
17
+ * The initial condition is u(t_0) = 1.0 and we use the parameters
18
+ * p = [1.5, 1.0, 3.0, 1.0]. The integration interval is t \in [0, 10.].
19
+ * The implicit BDF method from CVODES is used to solve the forward problem.
20
20
* Afterwards, the continuous adjoint sensitivity analysis capabilites of CVODES
21
- * are used to obtain the gradient of the cost function.
21
+ * are used to obtain the gradient of the cost function,
22
22
*
23
- * g(u,p,t ) = (sum(u)^2) / 2,
23
+ * g(u(t_f), p ) = (sum(u)^2) / 2
24
24
*
25
25
* with respect to the initial condition and the parameters.
26
26
* -----------------------------------------------------------------------------
@@ -89,8 +89,28 @@ void dgdu(N_Vector uvec, N_Vector dgvec)
89
89
dg [1 ] = u [0 ] + u [1 ];
90
90
}
91
91
92
+ // /* Function to compute the adjoint ODE right-hand side:
93
+ // -lambda*(df/du) - (dg/du)
94
+ // */
95
+ // static int adjoint_rhs(sunrealtype t, N_Vector uvec, N_Vector lvec,
96
+ // N_Vector ldotvec, void* user_data)
97
+ // {
98
+ // sunrealtype* p = (sunrealtype*)user_data;
99
+ // sunrealtype* u = N_VGetArrayPointer(uvec);
100
+ // sunrealtype* lambda = N_VGetArrayPointer(lvec);
101
+ // sunrealtype* lambdaDot = N_VGetArrayPointer(ldotvec);
102
+
103
+ // N_Vector dg = N_VClone(uvec);
104
+ // vjp(lvec, ldotvec, t, uvec, user_data);
105
+ // dgdu(uvec, dg);
106
+ // N_VLinearSum(-1.0, ldotvec, -1.0, dg, ldotvec);
107
+ // N_VDestroy(dg);
108
+
109
+ // return 0;
110
+ // }
111
+
92
112
/* Function to compute the adjoint ODE right-hand side:
93
- -lambda *(df/du) - (dg /du)
113
+ -mu *(df/du)
94
114
*/
95
115
static int adjoint_rhs (sunrealtype t , N_Vector uvec , N_Vector lvec ,
96
116
N_Vector ldotvec , void * user_data )
@@ -100,26 +120,6 @@ static int adjoint_rhs(sunrealtype t, N_Vector uvec, N_Vector lvec,
100
120
sunrealtype * lambda = N_VGetArrayPointer (lvec );
101
121
sunrealtype * lambdaDot = N_VGetArrayPointer (ldotvec );
102
122
103
- N_Vector dg = N_VClone (uvec );
104
- vjp (lvec , ldotvec , t , uvec , user_data );
105
- dgdu (uvec , dg );
106
- N_VLinearSum (-1.0 , ldotvec , -1.0 , dg , ldotvec );
107
- N_VDestroy (dg );
108
-
109
- return 0 ;
110
- }
111
-
112
- /* Function to compute the adjoint ODE right-hand side:
113
- -lambda*(df/du)
114
- */
115
- static int adjoint_rhs_alt (sunrealtype t , N_Vector uvec , N_Vector lvec ,
116
- N_Vector ldotvec , void * user_data )
117
- {
118
- sunrealtype * p = (sunrealtype * )user_data ;
119
- sunrealtype * u = N_VGetArrayPointer (uvec );
120
- sunrealtype * lambda = N_VGetArrayPointer (lvec );
121
- sunrealtype * lambdaDot = N_VGetArrayPointer (ldotvec );
122
-
123
123
vjp (lvec , ldotvec , t , uvec , user_data );
124
124
N_VScale (-1.0 , ldotvec , ldotvec );
125
125
0 commit comments