Group Abstract

Message Boards

WOLFRAM COMMUNITY

6.7K Views

2 Replies

1 Total Like

View groups...

Follow this post

Share this post:

GROUPS:

Physics Mathematica Calculus Equation Solving Graphics and Visualization Wolfram Language

plotting motion of two charged particles

Hayashi Yoshiaki

Posted 12 years ago

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]

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]

POSTED BY: Hayashi Yoshiaki

2 Replies

Sort By:

Hayashi Yoshiaki

Posted 12 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]

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 BY: Hayashi Yoshiaki

Nasser M. Abbasi

Nasser M. Abbasi, student

Posted 12 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"}]

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"}]

POSTED BY: Nasser M. Abbasi

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Feedback