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\limits_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