# plotting motion of two charged particles

Posted 9 years ago
5673 Views
|
2 Replies
|
1 Total Likes
|
 I wanted to plot motions of two charged particle by solving equation of motion with lorentz force so I write the code below referencing this http://www.wolfram.com/mathematica/new-in-9/advanced-hybrid-and-differential-algebraic-equations/double-pendulum.html but this don't show anything. someone please help!!!Mathematica code deqns2 = { m1 x1''[t] == (q1 q2)/(4 Pi e0 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2 + (z1[t] - z2[t])^2)^(3/2)) (x1[t] - x2[t] + (y2'[t] (z1[t] - z2[t]) - z2'[t] (y1[t] - y2[t]))), m1 y1''[t] == (q1 q2)/(4 Pi e0 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2 + (z1[t] - z2[t])^2)^(3/2)) (y1[t] - y2[t] + (z2'[t] (x1[t] - x2[t]) - x2'[t] (z1[t] - z2[t]))), m1 z1''[t] == (q1 q2)/(4 Pi e0 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2 + (z1[t] - z2[t])^2)^(3/2)) (z1[t] - z2[t] + (x2'[t] (y1[t] - y2[t]) - y2'[t] (x1[t] - x2[t]))), m2 x2''[t] == (q1 q2)/(4 Pi e0 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2 + (z1[t] - z2[t])^2)^(3/2)) (x2[t] - x1[t] + (y1'[t] (z2[t] - z1[t]) - z1'[t] (y2[t] - y1[t]))), m2 y2''[t] == (q1 q2)/(4 Pi e0 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2 + (z1[t] - z2[t])^2)^(3/2)) (y2[t] - y1[t] + (z1'[t] (x2[t] - x1[t]) - x1'[t] (z2[t] - z1[t]))), m2 z2''[t] == (q1 q2)/(4 Pi e0 ((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2 + (z1[t] - z2[t])^2)^(3/2)) (z2[t] - z1[t] + (x1'[t] (y2[t] - y1[t]) - y1'[t] (x2[t] - x1[t])))} ics2 = {x1[0] == 0, y1[0] == 0, z1[0] == 0, x1'[0] == 1, y1'[0] == 0, z1'[0] == 0,x2[0] == 1, y2[0] == 1, z2[0] == 1, x2'[0] == 0, y2'[0] == 1, z2'[0] == 0} params2 = {q1 -> 1, q2 -> 1, m1 -> 0.1, m2 -> 0.1, e0 -> 1} soldp2 = First[NDSolve[{deqns2, ics2} /. params2,{x1, y1, z1, x2, y2, z2}, {t, 0, 500}]] ParametricPlot3D[Evaluate[{x1[t], y1[t], z1[t]} /. soldp],{t, 0, 500}, ImageSize -> Medium] 
2 Replies
Sort By:
Posted 9 years ago
 Thank you myself! I solved my problem by myself. below is the thing I wanted to do. wierd motion of 2 charged particles!Mathematica Code r1 = {x1[t], y1[t], z1[t]}; r2 = {x2[t], y2[t], z2[t]}; v1 = Dt[r1, t]; v2 = Dt[r2, t]; deqns = Flatten[ { Thread[ m1 Dt[v1, t] == (q1 q2)/( 4 Pi e0 ) ((r1 - r2)/(Norm[r1 - r2])^3 + 1/c^2 (v1\[Cross](v2\[Cross]((r1 - r2)/(Norm[r1 - r2])^3))))], Thread[ m2 Dt[v2, t] == (q1 q2)/( 4 Pi e0 ) ((r2 - r1)/(Norm[r1 - r2])^3 + 1/c^2 (v2\[Cross](v1\[Cross](( r2 - r1)/(Norm[r1 - r2])^3))))] } ]; ics = {x1[0] == 0, y1[0] == 0, z1[0] == 0, x1'[0] == 0, y1'[0] == 0, z1'[0] == 0, x2[0] == 5, y2[0] == 0, z2[0] == 0, x2'[0] == 0.1, y2'[0] == 0.1, z2'[0] == 0.1}; params = {q1 -> 10, q2 -> -10, m1 -> 1, m2 -> 1, e0 -> 1 , c -> 1}; soldp = First[NDSolve[{deqns, ics} /. params, {x1, y1, z1, x2, y2, z2}, {t, 0, 50}, MaxSteps -> Infinity]]; ParametricPlot3D[ Evaluate[{{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}} /. soldp], {t, 0, 50}, PlotStyle -> {Red, Blue}, ImageSize -> Medium, PlotLegends -> {"Trajectory1", "Trajectory2"}] Animate[Graphics3D[{{PointSize[.025], {Red, Point[{x1[t], y1[t], z1[t]}]}, {Blue, Point[{x2[t], y2[t], z2[t]}]}} /. soldp, {Red, Line[ Map[Function[Evaluate[{x1[#], y1[#], z1[#]} /. soldp]], Range[0, t, .025]]]}, {Blue, Line[Map[Function[Evaluate[{x2[#], y2[#], z2[#]} /. soldp]], Range[0, t, .025]]]} }, PlotRange -> {{0, 10}, {0, 10}, {0, 10}}, Axes -> True, Ticks -> False, ImageSize -> 500], {t, 0, 50, .025}, SaveDefinitions -> True] 
Posted 9 years ago
 You did not copy it right. It works for me: deqns = {Subscript[m, 1] x1''[t] == (\[Lambda]1[t]/Subscript[l, 1]) x1[t] - (\[Lambda]2[t]/Subscript[l, 2]) (x2[t] - x1[t]), Subscript[m, 1] y1''[t] == (\[Lambda]1[t]/Subscript[l, 1]) y1[t] - (\[Lambda]2[t]/Subscript[l, 2]) (y2[t] - y1[t]) - Subscript[m, 1] g, Subscript[m, 2] x2''[t] == (\[Lambda]2[t]/Subscript[l, 2]) (x2[t] - x1[t]), Subscript[m, 2] y2''[t] == (\[Lambda]2[t]/Subscript[l, 2]) (y2[t] - y1[t]) - Subscript[m, 2] g}; aeqns = {x1[t]^2 + y1[t]^2 == Subscript[l, 1]^2, (x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 == Subscript[l, 2]^2}; ics = {x1[0] == 1, y1[0] == 0, x1'[0] == 0, y1'[0] == 0, x2[0] == 1, y2[0] == -1, x2'[0] == 0, y2'[0] == 0}; params = {g -> 9.81, Subscript[m, 1] -> 1, Subscript[m, 2] -> 1, Subscript[l, 1] -> 1, Subscript[l, 2] -> 1}; soldp = First[NDSolve[{deqns, aeqns, ics} /. params, {x1, y1, x2, y2, \[Lambda]1, \[Lambda]2}, {t, 0, 15}, Method -> {"IndexReduction" -> {Automatic, "IndexGoal" -> 0}}]]; ParametricPlot[Evaluate[{{x1[t], y1[t]}, {x2[t], y2[t]}} /. soldp], {t, 0, 15}, PlotStyle -> {Red, Blue}, ImageSize -> Medium, PlotLegends -> {"Trajectory of pendulum 1", "Trajectory of pendulum 2"}] 
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments