Hi there,

I've got a question to ask. 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) = \integrate\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)

Partial Derivatives:-

delta H/delta u = δH/δu = 2u+p2 =0

u=-p2/2;

- delta H/delta x1 = - δH/δx1 = -2x1;

- delta H/delta x2 = - δ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]; // all the derivatives

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 (newt) without using shooting method. Can anyone help me?

Thank you.

suli