Message Boards

WOLFRAM COMMUNITY

12782 Views

5 Replies

1 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Mathematica

Hamiltonian System in Mathematica

Phil R

Posted 11 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}]

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

Sort By:

Phil R

Posted 11 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

Ivan Morozov

Ivan Morozov, Budker Institute of Nuclear Physics

Posted 11 years ago

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

Paritosh Mokhasi

Paritosh Mokhasi, Wolfram Research Inc.

Posted 11 years ago

Phi, What are the definitions of Uff and coordreticolo?

POSTED BY: Paritosh Mokhasi

Phil R

Posted 11 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.109382610^(-31); q = -1.6021765310^(-19); Uff[r_] := ((qQ)/(4Piepsilon)) 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]}];

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

Paritosh Mokhasi

Paritosh Mokhasi, Wolfram Research Inc.

Posted 11 years ago

Looking at the first three equations in your system, they are defined as: Derivative[1][x][t] == -3.2933110^30 pxx[t] + Derivative[1][x][t], Derivative[1][y][t] == -3.2933110^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

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback