0
|
11224 Views
|
5 Replies
|
1 Total Likes
View groups...
Share
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}] 
5 Replies
Sort By:
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 10 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 10 years ago
 Phi, What are the definitions of Uff and coordreticolo?
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 10 years ago
 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.