Hamiltonian Question

[app08sfs]

Hi there,

I've got a question to ask regarding Hamiltonian problem. I'm doing on a simple question below:

State equations ODE:-

dy/dx (t) = x2(t)

d^2y/dx^2 (t) = -x2(t) + u(t) ....where u is a control.

Initial condition x1(0)=0 and x2(0)=0.


We try to minimize J(u) = \int_0^10... [y1^2(t) + u^2(t)] dt;

At the final time T=10, x1(10) and x2(10) is free.


Therefore the Hamiltonian is:-

H= x1^2 + u^2 + p1*x2 + p2*(-x2 + u)

δH/δu= 2u+p2 =0
u=-p2/2;

- δH/δx1= -2x1;

- δH/δx2= -(p1-p2);

I want p1(10) and p2(10) to be zero at final time T=10.

I've manage to solve it with shooting method (shoot) along with newton method (newt) below:

#include "stdio.h"
#include "nr3.h"
#include "ludcmp.h"
#include "stepper.h"
#include "odeint.h"
#include "qrdcmp.h"
#include "roots_multidim.h"
#include "stepperdopr5.h"
#include "stepperdopr853.h"
#include "shoot.h"

struct Rhs { // odes and partial derivative
Int m;
Doub c2;
Rhs() {}
void operator() (const Doub t, VecDoub_I &y, VecDoub_O &dydx)
{
Doub u= -y[3]/2.0;

dydx[0]= y[1];
dydx[1]= -y[1] + u;
dydx[2]= -2*y[0]; //p1=y[2]
dydx[3]= -(y[2]-y[3]); //p2=y[3]
}
};

struct Load { // initial conditions at t = x1 initial time
VecDoub y;
Load() : y(4) {}
VecDoub operator() (const Doub x1, VecDoub_I &v)
{
y[0]= 1.00;
y[1]= 0.00;
y[2]= v[0]; //guessing the value for y[2]
y[3]= v[1]; //guessing the value for y[3]

return y;
}
};

struct Score { // desired value(s) at t = x2 the endpoint
VecDoub f;
Score() : f(2) {}
VecDoub operator() (const Doub x2, VecDoub_I &y)
{
f[0]= y[2];
f[1]= y[3];

cout << "score = "<< f[0] << " " << f[1] << " y1[T] = " << y[0] << " y2[T] = " << y[1] << " p1[T] = " << y[2] << " p2[T] = " << y[3] << endl;
return f;
}
};


Int main() {
const Int N2=2; // number of guessed values at t = x1 //,MM=1;
Bool check;
VecDoub v(N2);
Int nvar=4; // number of variables

v[0]= 10.00; // initial guess value for y[2]
v[1]= 10.00; // guess value for y[3]

Doub x1=0.0; // initial time
Doub x2=1.0; // final time

Load load; // initial values

Rhs rhs; // odes

Score score; // desired endpoint values

Shoot<Load,Rhs,Score> shoot(nvar,x1,x2,load,rhs,score);
newt(v,check,shoot);

system("PAUSE");
return 0;
}


Now the problem is that I want to solve the same problem just with Newton Method without using shooting method. Can anyone help me?

Thank you.

suli