# Cannot solve differential equation for electron motion in Electric Field

Posted 9 years ago
4137 Views
|
|
3 Total Likes
|
 Hi everyone, I'm trying to find the trajectory of an electron in a graphene lattice. I've a set of three differential equations in the three directions mx''=Fx,my''=Fy,mz''=Fz and i try to solve my system whith NDSolve and to plot the result, but what I have in output is the same linear function whatever are boundary conditions. It's clear that it's not a right solution, so could you help me to understand what I'm doing wrong? Here's my code: Clear["Global*"] coordx = {0.6327, 1.88058, 3.03927, 4.28716, 5.44584, 6.69373, 7.85241, 9.10029, 1.9728, 3.22069, 4.37937, 5.62726, 6.78594, 8.03382, 9.19251, 10.4404, 3.3129, 4.56079, 5.71947, 6.96736, 8.12604, 9.37393, 10.53261, 11.7805, 4.653, 5.90089, 7.05956, 8.30746, 9.46614, 10.71403, 11.87271, 13.1206}; coordy = {1.31579, 1.79441, 0.98308, 1.46168, 0.65035, 1.12896, 0.31763, 0.79624, 3.08485, 3.56345, 2.75213, 3.23074, 2.41941, 2.89802, 2.08669, 2.56529, 4.8539, 5.33251, 4.52119, 4.99979, 4.18846, 4.66707, 3.85574, 4.33435, 6.62295, 7.10157, 6.29024, 6.76884, 5.95751, 6.43613, 5.62479, 6.10341}; coordz = {3.01288, 2.54181, 2.65281, 2.18174, 2.29275, 1.82168, 1.93268, 1.46161, 1.95888, 1.4878, 1.5988, 1.12774, 1.23874, 0.76767, 0.87868, 0.4076, 0.90487, 0.43379, 0.5448, 0.07373, 0.18473, -0.28634, -0.17533, -0.64641, -0.14914, -0.62022, \ -0.5092, -0.98028, -0.86927, -1.34035, -1.22934, -1.70042}; me = 9.01*10^-31; pi = 3.14159; epsilon0 = 8.854*10^-12; q = -1.6*10^-19; Q = 1.6*10^-19; step = 0.01; Forzax[p_, r_] := Sum[(Q*q)/(4 pi*epsilon0*Norm[r - p[[i]]]^2), {i, Length[p]}] Forzay[p_, r_] := Sum[(Q*q)/(4 pi*epsilon0*Norm[r - p[[i]]]^2), {i, Length[p]}] Forzaz[p_, r_] := Sum[(Q*q)/(4 pi*epsilon0*Norm[r - p[[i]]]^2), {i, Length[p]}] step = 0.01; sol = NDSolve[{x''[t] - Forzax[coordx, {x[t]}]/me == 0, y''[t] - Forzay[coordy, {y[t]}]/me == 0, z''[t] - Forzaz[coordz, {z[t]}]/me == 0, x[0] == 0, y[0] == -5, z[0] == -4, x'[0] == 1, y'[0] == 2, z'[0] == 0}, {x[t], y[t], z[t]}, {t, 0, 15000}] Show[ParametricPlot3D[ Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, 15000}, PlotRange -> All]] `Thank you in advance. Regards Filippo