Message Boards Message Boards

0
|
11587 Views
|
5 Replies
|
1 Total Likes
View groups...
Share
Share this post:
GROUPS:

Hamiltonian System in Mathematica

Posted 10 years ago

Hi everyone, I'm trying to solve the Hamiltonian system to find the trajectory of an electron in a lattice. I've found the equations of my system but i can't solve the system. NDSOlve give me that error:

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations. >>
NDSolve::index: The DAE solver failed at t = 0.`. The solver is intended for index 1 DAE systems and structural analysis indicates that the DAE index is 2. The option Method->{"IndexReduction"->Automatic} may be used to reduce the index of the system. >>

Trying with Method->{"IndexReduction"->Automatic} option it gives me:

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations. >>
NDSolve::ndcf: Repeated convergence test failure at t == 0.`; unable to continue. >>

(I've not used the LegandreTransform function because this work is for academic purpose) Could someone help me please? Thank you in advance, here is my code:

L = 1/(2 m) (pxx[t]^2 + pyy[t]^2 + pzz[t]^2) - Uff[coordreticolo];
H[p_, q_, dq_] := Sum[p[[i]]*dq[[i]] - L, {i, Length[q]}];
H[{pxx[t], pyy[t], pzz[t]}, {x[t], y[t], z[t]}, {x'[t], y'[t], z'[t]}]
xp = D[H[{pxx[t], pyy[t], pzz[t]}, {x[t], y[t], z[t]}, {x'[t], y'[t], 
    z'[t]}], pxx[t]]
yp = D[H[{pxx[t], pyy[t], pzz[t]}, {x[t], y[t], z[t]}, {x'[t], y'[t], 
    z'[t]}], pyy[t]]
zp = D[H[{pxx[t], pyy[t], pzz[t]}, {x[t], y[t], z[t]}, {x'[t], y'[t], 
    z'[t]}], pzz[t]]
pxp = -D[H[{pxx[t], pyy[t], pzz[t]}, {x[t], y[t], z[t]}, {x'[t], 
     y'[t], z'[t]}], x[t]]
pyp = -D[H[{pxx[t], pyy[t], pzz[t]}, {x[t], y[t], z[t]}, {x'[t], 
     y'[t], z'[t]}], y[t]]
pzp = -D[H[{pxx[t], pyy[t], pzz[t]}, {x[t], y[t], z[t]}, {x'[t], 
     y'[t], z'[t]}], z[t]]
NDSolve[{x'[t] == xp, y'[t] == yp, z'[t] == zp, pxx'[t] == pxp, 
  pyy'[t] == pyp, pzz'[t] == pzp, x[0] == 0, pxx[0] == 0, y[0] == 0, 
  pyy[0] == 0, z[0] == 0, pzz[0] == 0}, {x[t], y[t], z[t], pxx[t], 
  pyy[t], pzz[t]}, {t, 0, 15}, 
 Method -> {"IndexReduction" -> Automatic}]
POSTED BY: Phil R
5 Replies
Posted 10 years ago

Thank you Paritosh and Ivan. I have expressed L through p because of H(q,p,t)=Sum[q'[t]*p[t]] -L[q,p,t] so I wanted to solve my system in x,y,z,px,py,pz. "IndexReduction" should not be used, Mathematica gave me that suggestion because my code was uncorrect. With the changes that you suggested, and having reinitialized the variables x'[t],y'[t] and z'[t] now it seems to work.

Regards Filippo

POSTED BY: Phil R

Hi, Phil,

This line is wrong:

  H[p_, q_, dq_] := Sum[p[[i]]*dq[[i]] - L, {i, Length[q]}];

should be:

H[p_, q_, dq_] := Sum[p[[i]]*dq[[i]] , {i, Length[q]}] - L;

Also, why $L$ is expressed through $p_i$? Should "IndexReduction" method be necessarily used? As it is seen from documentation "IndexReduction" is used mostly for constrained problems. You system is unconstrained and it is straightforward to obtain Hamiltonian.

I.M.

POSTED BY: Ivan Morozov

Phi, What are the definitions of Uff and coordreticolo?

POSTED BY: Paritosh Mokhasi
Posted 10 years ago

Sorry, you're right, I haven't posted that part of code for simplicity. Uff is lattice's ions potential; coordreticolo is matrix of ions positions(x,y,z). Here's the definition code:

coordreticolo = {{0.6327, 1.31579, 3.01288},
   {1.88058, 1.79441, 2.54181},
   {3.03927, 0.98308, 2.65281},
   {4.28716, 1.46168, 2.18174},
   {5.44584, 0.65035, 2.29275},
   {6.69373, 1.12896, 1.82168},
   {7.85241, 0.31763, 1.93268},
   {9.10029, 0.79624, 1.46161},
   {1.9728, 3.08485, 1.95888},
   {3.22069, 3.56345, 1.4878},
   {4.37937, 2.75213, 1.5988},
   {5.62726, 3.23074, 1.12774},
   {6.78594, 2.41941, 1.23874},
   {8.03382, 2.89802, 0.76767},
   {9.19251, 2.08669, 0.87868},
   {10.4404, 2.56529, 0.4076},
   {3.3129, 4.8539, 0.90487},
   {4.56079, 5.33251, 0.43379},
   {5.71947, 4.52119, 0.5448},
   {6.96736, 4.99979, 0.07373},
   {8.12604, 4.18846, 0.18473},
   {9.37393, 4.66707, -0.28634},
   {10.53261, 3.85574, -0.17533},
   {11.7805, 4.33435, -0.64641},
   {4.653, 6.62295, -0.14914},
   {5.90089, 7.10157, -0.62022},
   {7.05956, 6.29024, -0.5092},
   {8.30746, 6.76884, -0.98028},
   {9.46614, 5.95751, -0.86927},
   {10.71403, 6.43613, -1.34035},
   {11.87271, 5.62479, -1.22934},
   {13.1206, 6.10341, -1.70042}};
Q = 1.60217653 10^(-19);
epsilon = 8.85418781762 10^(-12);
m = 9.1093826*10^(-31);
q = -1.60217653*10^(-19);
Uff[r_] := ((q*Q)/(4*Pi*epsilon))*
   Sum[1/((x[t] - r[[i, 1]])^2 + (y[t] - r[[i, 2]])^2 + (z[t] - 
           r[[i, 3]])^2)^(1/2), {i, Length[r]}];
POSTED BY: Phil R

Looking at the first three equations in your system, they are defined as:

Derivative[1][x][t] == -3.29331*10^30 pxx[t] + Derivative[1][x][t], 
Derivative[1][y][t] == -3.29331*10^30 pyy[t] + Derivative[1][y][t], 
Derivative[1][z][t] == -3.29331*10^30 pzz[t] + Derivative[1][z][t]

This means that pxx, pyy, ptt are zero and so are their derivatives. This means that you end up with three algebraic equations for x,y,z from the last 3 equations. Since this is not my field of expertise, perhaps you can comment on whether the formulation of the equations are correct or not. If they are, then this is an algebraic problem and is better suited to be solved using FindRoot, NSolve.

POSTED BY: Paritosh Mokhasi
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract