Numerical solutions of differential equations of motion in three dimensions

Numerical solutions of differential equations of motion in three dimensions

The motion of a particle can be determined from the force acting on the particle. The relation between the force $\vec{F}$ and the position of the particle as a function of time $\vec{r}(t)$ is given by Newton's law,

$\large \vec{F}(\vec{r},\vec{v},t) = m\frac{d^2\vec{r}}{dt^2}$.

Generally, the force $\vec{F}$ can depend on the position $\vec{r}$, the velocity $\vec{v}$ and the time $t$. Here $m$ is the mass of the particle. Newton's law can be rewritten as six first order differential equations,

$\large \frac{dx}{dt}=v_x$, $\large \frac{dv_x}{dt}=F_x(x,y,z,v_x,v_y,v_z,t)/m$,

$\large \frac{dy}{dt}=v_y$, $\large \frac{dv_y}{dt}=F_y(x,y,z,v_x,v_y,v_z,t)/m$,

$\large \frac{dz}{dt}=v_z$, $\large \frac{dv_z}{dt}=F_z(x,y,z,v_x,v_y,v_z,t)/m$.

These equations can be solved numerically step-by-step. If the initial conditions for $\vec{x}(t_0)$ and $\vec{v}(t_0)$ are known, a good estimate for $\vec{x}$ and $\vec{v}$ a short time $\Delta t$ later is,

$\large \vec{x}(t_0+\Delta t) \approx \frac{d\vec{v}}{dt}|_{t_0}\Delta t$ and $\large \vec{v}(t_0+\Delta t) \approx \vec{F}(x(t_0),v_x(t_0),t_0)\Delta t/m.$

Once an estimate for the position and the velocity of the object at time $t_0 + \Delta t$ is calculated, they can be used to estimate the position and the velocity at time $t_0 + 2\Delta t$. By repeating this over and over, a table of times $(t)$, positions $(x,y,z)$, velocities $(v_x,v_y,v_z)$, and forces $(F_x,F_y,F_z)$ can be made. Other quantities such as the kinetic energy, the work performed, and the instantaneous power can be calculated. The kinetic energy is $E_{\text{kin}}=\frac{v_x^2+v_y^2+v_z^2}{2m}$ [J]. When a force moves a particle a distance $d\vec{r}$ it performs an increment of work $dW=\int\vec{F}\cdot d\vec{r}$ [J]. This work is performed in a time $dt$. The power needed for this is $P=dW/dt=\vec{F}\cdot\vec{v}$ [W].

In the form below, $F_x$, $F_y$, and $F_z$ are the three components of the force, $m$ is the mass, and $t$ is the time. Starting from intitial conditions, the equations are integrated for a total number of $N_{steps}$ steps using a step size of $\Delta t$. The force can be specified as a function of the positions $x$, $y$, $z$, the velocities $v_x$, $v_y$, $v_z$ and the time $t$.

= , =

the animation to zoom or rotate.

$t$ [s] $x$ [m] $y$ [m] $z$ [m] $v_x$ [m/s] $v_y$ [m/s] $v_z$ [m/s] $F_x$ [N] $F_y$ [N] $F_z$ [N] $P$ [W] $E_{\text{kin}}$ [J] $W$ [J]

Questions

<p>An object moving in three dimensions is described by six variables: $x$, $y$, $z$, $v_x$, $v_y$, and $v_z$. If the force on the object is known, then the motion can be described by six first order differential equations,</p>

<p style="text-align:center">$\large \frac{dx}{dt}=v_x$   $\large \frac{dv_x}{dt}=F_x(x,y,z,v_x,v_y,v_z,t)/m$</p>
<p style="text-align:center">$\large \frac{dy}{dt}=v_y$   $\large \frac{dv_y}{dt}=F_y(x,y,z,v_x,v_y,v_z,t)/m$</p>
<p style="text-align:center">$\large \frac{dz}{dt}=v_z$   $\large \frac{dv_z}{dt}=F_z(x,y,z,v_x,v_y,v_z,t)/m$</p>

<p>These equations can be solved numerically step-by-step. If the initial conditions for $\vec{x}(t_0)$ and $\vec{v}(t_0)$ are known, a good estimate for $\vec{x}$ and $\vec{v}$ a short time $\Delta t$ later is,</p>

<p style="text-align:center">$\large \vec{x}(t_0+\Delta t) \approx \frac{d\vec{v}}{dt}|_{t_0}\Delta t$ and $\large \vec{v}(t_0+\Delta t) \approx \vec{F}(x(t_0),v_x(t_0),t_0)\Delta t/m.$</p>

<p>Once an estimate for the position and the velocity of the object at time $t_0 + \Delta t$ is calculated, they can be used to estimate the position and the velocity at time $t_0 + 2\Delta t$. By repeating this over and over, a table of times $(t)$, positions $(x,y,z)$, velocities $(v_x,v_y,v_z)$, and forces  $(F_x,F_y,F_z)$ can be made. Other quantities such as the kinetic energy, the work performed, and the instantaneous power can be calculated. The kinetic energy is $E_{\text{kin}}=\frac{v_x^2+v_y^2+v_z^2}{2m}$ [J]. When a force moves a particle a distance $d\vec{r}$ it performs an increment of work $dW=\int\vec{F}\cdot d\vec{r}$ [J]. This work is performed in a time $dt$. The power needed for this is $P=dW/dt=\vec{F}\cdot\vec{v}$ [W]. </p>

<p>In the form below, $F_x$, $F_y$, and $F_z$ are the three components of the force, $m$ is the mass, and $t$ is the time. Starting from intitial conditions, the equations are integrated for a total number of $N_{steps}$ steps using a step size of $\Delta t$. The force can be specified as a function of the positions $x$, $y$, $z$, the velocities $v_x$, $v_y$, $v_z$ and the time $t$.</p>

<h3 style="text-align:center">Satellite orbits</h3>

<p>A satellite orbits the earth. The graviational force on the satellite is,</p>

<p style="text-align:center">$\large \vec{F} = -\frac{Gm_e m_s}{r^2} \hat{r},$</p>

<p>where $G= 6.6726 \times 10^{-11}$ N m²/kg² is the graviational constant, $m_e = 5.97219 \times 10^{24}$ kg is the mass of the earth, $m_s$ is the mass of the satellite, and $\vec{r}$ is the position of the satellite measured from the center of the earth.</p>

<form name="sat">
<p>As long as the earth is much heavier than the satellite, the mass of the satellite does not matter. The satellite mass does not appear in the expression for the acceleration. Only the intitial position and velocity of the satellite determines its orbit. If the initial velocity is above the escape velocity, the satellite will escape from the earth's gravitational field and travel away across the galaxy.</p>
<p>The initial conditions at $t=0$ are: <br />
$x=$ <input type="text" size="10" value="0" name="x"> m  $y=$ <input type="text" size="10" value="6371000" name="y"> m  $z=$ <input type="text" size="10" value="0" name="z"> m 
$v_x=$ <input type="text" size="10" value="7900" name="vx"> m/s  $v_y=$ <input type="text" size="10" value="0" name="vy"> m/s  $v_z=$ <input type="text" size="10" value="0" name="vz"> m/s  $m_s=$ <input type="text" size="10" value="100" name="mass"> kg</p>

<p style="text-align:center">
<input type="button" value="Load these equations into the form below" onclick='satellite();'>
</p>
</form>

<p>If the orbit falls below about 6400000 m, the satellite will crash into the earth. There are various kinds of orbits such as geosynchronous orbits, geostationary orbits, low earth orbits, elliptical orbits, and graveyard orbits. The difference just depends on the initial conditions of the satellite. A long time step should be used to calculate satelite orbits. </p>

<h3 style="text-align:center">Inclined plane</h3>

<p>A block of mass $m$ is put on an inclined plane that makes an angle of $\theta$ degrees from horizontal. </p>

<table align="center"><tr><td>$y$</td><td><div id="placeholder_block" style="width:300px;height:300px;"></div></td></tr>
<tr><td></td><td><p style="text-align:center">$x$</p></td></tr></table>

<form name="ball">

<table align="center">

<tr><td><p style="text-align:right; margin-right:0cm;">$\theta=$ <span id="btext">45</span> °</td><td><input type="button" value="-" onclick="ball.slide_b.value=ball.slide_b.value-1; inclined();"><input type='range' min='0' max='90' step='1' name='slide_b' value='45' onchange='inclined();' oninput='inclined();'/><input type="button" value="+" onclick="ball.slide_b.value=eval(ball.slide_b.value) + 1; inclined ();"></td></tr>

<tr><td><p style="text-align:right; margin-right:0cm;">$m=$ <span id="mtext">0.4</span> [kg]</td><td><input type="button" value="-" onclick="ball.slide_m.value=ball.slide_m.value-0.01; inclined();"><input type='range' min='0.01' max='1' step='0.01' name='slide_m' value='0.4' onchange='inclined();' oninput='inclined();'/><input type="button" value="+" onclick="ball.slide_m.value=eval(ball.slide_m.value) + 0.01; inclined();"></td></tr>

</table>

</form>

<p>The gravitational force, $-mg\hat{y}$ can be decomposed into a force normal to the inclined plane and a force parallel to the inclined plane.</p>

<p style='text-align:center'><img src='block1.jpg'></p>

<p>The force parallel to the plane, $|F_{\parallel}|=mg\sin\theta$, causes the block to slide along the plane. In terms of its $x$- and $y$-components, the force along the plane is,</p>

<p style='text-align:center'><img src='block2.jpg'></p>

<p style='text-align:center'>$\vec{F}_{\parallel} = -mg\sin\theta\cos\theta \hat{x} - mg\sin\theta\sin\theta \hat{y}.$</p>

<p>This is a constant force; it is independent of the position, the velocity, and the time. We expect the mass to exhibit parabolic motion and formulas for the position vector and the velocity vector can be found on the page <a href="http://lampz.tugraz.at/~hadley/physikm/apps/rvaf/constant_force.en.php">Constant Force = Parabolic Motion</a>. It is of course also possible to put this force and the initial conditions in a numerical differential equation solver. The force and the mass is loaded into the differential solver below every time the sliders are adjusted. A drag force can be added to the differential equation solver by adding a term to the force that is in the opposite direction of the velocity. If this drag force is linear $\vec{F}_{drag} = -a\vec{v}$, then it is possible to find <a href="http://lampz.tugraz.at/~hadley/physikm/apps/rvaf/terminal.en.php">formulas to describe the resulting motion</a>. However, if the drag force is nonlinear such as, $\vec{F}_{drag} = -a|\vec{v}|\vec{v}$, no simple formulas can be found and a numerical differential equation solver is the best option.</p>

<p style="text-align:center">No drag force: $\vec{F} =  -mg\sin\theta\cos\theta \,\hat{x} - mg\sin\theta\sin\theta \,\hat{y}$<br />
Linear drag force: $\vec{F} =  -mg\sin\theta\cos\theta\, \hat{x} - mg\sin\theta\sin\theta \,\hat{y}-a\vec{v}$<br />
Nonlinear drag force: $\vec{F} =  -mg\sin\theta\cos\theta \,\hat{x} - mg\sin\theta\sin\theta \,\hat{y}-a|\vec{v}|\vec{v}$</p>

<p style="text-align:center"><input type="button" value="Add a linear drag force" onclick="linear();"> <input type="button" value="Add a nonlinear drag force" onclick="nonlinear();"></p>

<h3 style="text-align:center">Air friction</h3>

<p>A stone of mass 100 g is thrown with an intitial velocity $|\vec{v}|$ at an angle $\theta$ from horizontal. The stone is released from position $\vec{r}=0$. If air friction is ignored, the equation that describes the motion of the ball is, </p>

<p style="text-align:center">$\vec{r} = vt\cos\theta \,\hat{x}+ \left( vt\sin\theta -\frac{gt^2}{2}\right)\,\hat{y}\,\text{m},$</p>

<p>where $g = 9.81$ m/s² is the acceleration of gravity at the earth's surface and $t$ is measured in seconds. The stone follows a parabola and returns to the height from which it was thrown, $y=0$, at a distance $x$ from where it was thrown. For a given initial velocity, the angle that maximizes the distance $x$ is $\theta = \frac{\pi}{4}\,\text{rad} = $ 45°. </p>

<p>If a linear friction force is included in the problem $\vec{F}_{drag} = -b\vec{v}$, then the optimal angle changes. The form below calculates the trajectory of the stone for different angles, friction contants, and intitial velocities.</p>

<p style="text-align:center">$\vec{F} = m\frac{d^2\vec{r}}{dt^2} =  -b\vec{v}-mg \,\hat{z}$</p>

<form name="ball">
<p style="text-align:center">$b=$ <input type="text" size="10" value="1" name="bball"> N s/m</p>

<p style="text-align:center">The initial conditions at time $t=0$ are $\vec{r}=0$ <br />

$|\vec{v}| =$ <input type="text" size="10" value="100" name="vball"> m/s  $\theta=$ <input type="text" size="10" value="45" name="theta"> deg
</p>
<p style="text-align:center">
<input type="button" value="Load these values into the form below" onclick="airfriction();">
</p>
</form>

<h3 style="text-align:center">A ball is thrown in the wind</h3>

<p>A ball of mass $m$ is thrown and experiences a frictional drag force as it moves through a gas or a fluid. The forces acting on this ball are gravity $-mg\hat{z}$ and the drag force. If a wind is blowing, the drag force can be described by,</p>

<p style="text-align:center">$\vec{F}_{fric} = -a(\vec{v}-\vec{v}_{\text{wind}}) - b(\vec{v}-\vec{v}_{\text{wind}})|(\vec{v}-\vec{v}_{\text{wind}})|,$</p>

<p>where $a$ and $b$ are constants and $\vec{v}_{\text{wind}}$ is the velocity of the wind which can depend on position and time. For low <a href="http://en.wikipedia.org/wiki/Reynolds_number">Reynolds number</a>, the linear term $-a(\vec{v}-\vec{v}_{\text{wind}})$ usually dominates whereas for high Reynolds number, the quadratic term $- b(\vec{v}-\vec{v}_{\text{wind}})|(\vec{v}-\vec{v}_{\text{wind}})|$ dominates. </p>

<p style="text-align:center">$\vec{F}= m\frac{d^2\vec{r}}{dt^2} = -a(\vec{v}-\vec{v}_{\text{wind}}) - b(\vec{v}-\vec{v}_{\text{wind}})|(\vec{v}-\vec{v}_{\text{wind}})|-mg\,\hat{z}$</p>

<form name="ball">
<p>$m=$ <input type="text" size="10" value="0.1" name="mball"> kg $a=$ <input type="text" size="10" value="0.01" name="aball"> N s/m $b=$ <input type="text" size="10" value="0.03" name="bball"> N s²/m²</p>
<p>The three components of the wind vectors can be functions of space and time.<br />
$v_{\text{wind},x}=$ <input style="border: solid 1px blue" type="text" size="10" value="1" name="vwx"> m/s $v_{\text{wind},y}=$ <input style="border: solid 1px blue" type="text" size="10" value="7*exp(-x*x)" name="vwy"> m/s $v_{\text{wind},z}=$ <input style="border: solid 1px blue" type="text" size="10" value="-3*exp(-t*t)" name="vwz"> m/s</p>
<p>The initial conditions at time $t=0$ are, <br />
$x=$ <input type="text" size="10" value="0" name="xball"> m  $y=$ <input type="text" size="10" value="0" name="yball"> m  $z=$ <input type="text" size="10" value="0" name="zball"> m 
$v_x=$ <input type="text" size="10" value="-7" name="vxball"> m/s  $v_y=$ <input type="text" size="10" value="5" name="vyball"> m/s  $v_z=$ <input type="text" size="10" value="10" name="vzball"> m/s<br />
</p>
<p style="text-align:center">
<input type="button" value="Load these equations into the form below" onclick="wind();">
</p>
</form>

<h3 style="text-align:center">Rocket launch</h3>

<p>A model rocket with a mass $m$ is launched with a motor that provides a upward force of $F_{\text{thrust}}$ for 3 seconds. This force can be expressed mathematically as $\vec{F}=F_{\text{thrust}}H(3-t)\hat{z}$ where $H(x)$ is the <a href="https://en.wikipedia.org/wiki/Heaviside_step_function">Heaviside step function</a>. As the rocket uses fuel, its mass decreases. In the example below, the rocket looses half of its mass in the three seconds that it accelerates and then the mass remains constant, $m=0.1(2-H(3-t)t/3-H(t-3))$.  Other forces that act on the rocket are gravity $-mg\hat{z}$ and a drag force that can be described by,</p>

<p style="text-align:center">$\large \vec{F}_{fric} = -a(\vec{v}-\vec{v}_{\text{wind}}) - b(\vec{v}-\vec{v}_{\text{wind}})|(\vec{v}-\vec{v}_{\text{wind}})|,$</p>

<p>where $a$ and $b$ are constants and $\vec{v}_{\text{wind}}$ is the velocity of the wind which can depend on position and time. </p>

<p style="text-align:center">$\vec{F}=0.1(2-H(3-t)t/3-H(t-3))\frac{d^2\vec{r}}{dt^2} = F_{\text{thrust}}H(3-t)\hat{z}-a(\vec{v}-\vec{v}_{\text{wind}}) - b(\vec{v}-\vec{v}_{\text{wind}})|(\vec{v}-\vec{v}_{\text{wind}})|-mg\,\hat{z}$</p>

<form name="ball">
<p>$m=$ <input type="text" size="10" value="0.1*(2-H(3-t)*t/3-H(t-3))" name="mball"> kg $F_{\text{thrust}}=$ <input type="text" size="10" value="5" name="thrust"> N $a=$ <input type="text" size="10" value="0.01" name="aball"> N s/m $b=$ <input type="text" size="10" value="0.03" name="bball"> N s²/m²</p>
<p>The three components of the wind vectors can be functions of space and time.<br />
$v_{\text{wind},x}=$ <input style="border: solid 1px blue" type="text" size="20" value="0.2*exp(-t*t/50)" name="vwx"> m/s $v_{\text{wind},y}=$ <input style="border: solid 1px blue" type="text" size="20" value="-0.3*exp(-t*t/50)" name="vwy"> m/s $v_{\text{wind},z}=$ <input style="border: solid 1px blue" type="text" size="20" value="0" name="vwz"> m/s</p>
<p>The initial conditions at time $t=0$ are, <br />
$x=$ <input type="text" size="10" value="0" name="xball"> m  $y=$ <input type="text" size="10" value="0" name="yball"> m  $z=$ <input type="text" size="10" value="0" name="zball"> m 
$v_x=$ <input type="text" size="10" value="0" name="vxball"> m/s  $v_y=$ <input type="text" size="10" value="0" name="vyball"> m/s  $v_z=$ <input type="text" size="10" value="0" name="vzball"> m/s<br />
</p>
<p style="text-align:center">
<input type="button" value="Load these equations into the form below" onclick="rocket();">
</p>
</form>

<h3 style="text-align:center">A bouncing ball</h3>

<p>A ball of mass $m$ is thrown and bounces on the floor. The forces acting on this ball are gravity $-mg\hat{z}$, a drag force proportional to the velocity$-a\vec{v}$, and the force exerted by the floor when it bounces. The floor can be described as an elastic material that pushes up with a force $-kzH(-z)\hat{z}$. Here $k$ is the elastic contant and the Heaviside function $H(-z)$ ensures that this force only if the ball falls below $z=0$. </p>

<p>The larger the elastic constant $k$ is, the less the ball will go below $z=0$ and the simulation will approach bounces off a hard floor. However, some care must be taken in choosing the time step $\Delta t$. It should be much smaller than $2\pi\sqrt{m/k}$ so that the solver can describe the bounces off the floor properly. Some differential equation solvers adjust the time step automatically to account for problems like this. </p>

<p style="text-align:center">$\vec{F} = m\frac{d^2\vec{r}}{dt^2} = -a\vec{v}-mg\hat{z} -kzH(-z)\hat{z}$</p>

<form name="ball">
<p>$m=$ <input type="text" size="10" value="0.1" name="mball"> kg $a=$ <input type="text" size="10" value="0.01" name="aball"> N s/m $k=$ <input type="text" size="10" value="100" name="bball"> N/m</p>

<p>The initial conditions at time $t=0$ are, <br />
$x=$ <input type="text" size="10" value="0" name="xball"> m  $y=$ <input type="text" size="10" value="0" name="yball"> m  $z=$ <input type="text" size="10" value="0" name="zball"> m 
$v_x=$ <input type="text" size="10" value="-7" name="vxball"> m/s  $v_y=$ <input type="text" size="10" value="5" name="vyball"> m/s  $v_z=$ <input type="text" size="10" value="10" name="vzball"> m/s<br />
</p>
<p style="text-align:center">
<input type="button" value="Load these equations into the form below" onclick="bounce();">
</p>
</form>

<h3 style="text-align:center">Curveball: Magnus effect</h3>

<table align="right"><tr><td> </td><td><img src="https://upload.wikimedia.org/wikipedia/commons/thumb/f/f0/Magnus-anim-canette.gif/330px-Magnus-anim-canette.gif"><br /><a href="https://commons.wikimedia.org/wiki/File:Magnus-anim-canette.gif">Wikimedia Commons</a><br /><br />
</td></tr></table>

<p>A spinning ball thrown through the air will follow a curved path due to the <a href="https://en.wikipedia.org/wiki/Magnus_effect">Magnus effect</a>. The Magnus force on a smooth ball is,</p>

<p style="text-align:center">$\vec{F}_{\text{Magnus}} = \frac{4}{3} \pi\rho r^3(\vec{\omega}\times\vec{v}),$</p>

<p>where $r$ is the radius of the ball, $\rho$ is the density of the air, $\vec{\omega}$ is the angular frequency of the rotation, and $\vec{v}$ is the velocity of the ball. The drag force in air is approximated as,</p>

<p style="text-align:center">$\vec{F}_{\text{drag}} = -\frac{\pi}{2} r^2 \rho C_d|\vec{v}|\vec{v},$</p>

<p>where $C_d$ is the <a href="https://www1.grc.nasa.gov/beginners-guide-to-aeronautics/drag-of-a-sphere/">drag coefficient</a>. </p>

<p style="text-align:center">$\vec{F}= m\frac{d^2\vec{r}}{dt^2} = -\frac{\pi}{2} r^2 \rho C_d|\vec{v}|\vec{v}+ \frac{4}{3} \pi\rho r^3(\vec{\omega}\times\vec{v})-mg\,\hat{z}$</p>

<form name="ball">
<p>$m=$ <input type="text" size="10" value="0.1" name="mball"> kg $r=$ <input type="text" size="10" value="0.1" name="rball"> m</p>

<p>$\rho=$ <input type="text" size="10" value="1.293" name="rho"> kg/m³ $C_d=$ <input  style="border: solid 1px blue" type="text" size="20" value="0.47" name="Cd"></p>

<p>The three components of the vector that describes the rotation:<br />
$\vec{\omega}=$ <input style="border: solid 1px blue" type="text" size="20" value="0" name="wx"> $\hat{x}$  + <input style="border: solid 1px blue" type="text" size="20" value="0" name="wy"> $\hat{y}$  + <input style="border: solid 1px blue" type="text" size="20" value="10*exp(-0.001*t)" name="wz"> $\hat{z}$ rad/s</p>

<p>Generally, the drag coefficient is a function of the velocity and can be expressed in terms of vx, vy, and vz. The rotational velocity will decay with time due to friction. It can be expressed as a function of time.</p>

<p>The initial conditions at time $t=0$ are, <br />
$x=$ <input type="text" size="10" value="0" name="xball"> m  $y=$ <input type="text" size="10" value="0" name="yball"> m  $z=$ <input type="text" size="10" value="0" name="zball"> m 
$v_x=$ <input type="text" size="10" value="30" name="vxball"> m/s  $v_y=$ <input type="text" size="10" value="0" name="vyball"> m/s  $v_z=$ <input type="text" size="10" value="10" name="vzball"> m/s<br />
</p>
<p style="text-align:center">
<input type="button" value="Load these equations into the form below" onclick="magnus();">
</p>
</form>

<h3 style="text-align:center">Electric deflection of electrons</h3>

<p>Electrons are accelerated through a voltage $V_x$ towards a positively charged plate. Some of the electrons pass through a small hole in the plate and form and electron beam that travels to a region where an electric field is established by applying a voltage $V_y$ between two metal plates spaced a distance $d$ apart.</p>

<p style="text-align:center"><img src="../electricity/electric_deflection.png"></p>

<p>The electrons are accelerated from rest with a voltage of $V_x=5000$ V to the positively charged plate in a distance of 5 cm. The electrons get deflected as they pass through the electric field between the plates where $V_y = 60$ V. The force in this case is,</p>

<p style="text-align:center">$\vec{F}=m\frac{d^2\vec{r}}{dt^2} = qE_xH(0.05-x)\,\hat{x} +qE_yH(0.1-|x-0.25|)\,\hat{y}$.</p>

<p>A similation of electric deflection is found at <a href="http://lampz.tugraz.at/~hadley/physikm/script/electricity/electric_deflection_F.en.php">Electric deflection of an electron beam</a>.</p>

<h3 style="text-align:center">A charged particle in a constant magnetic field</h3>

<p>The force on a particle with charge $q$ in an magnetic field $\vec{B}$ is,</p>

<p style="text-align:center">$\vec{F} =q\vec{v}\times\vec{B}.$</p>

<p>If the velocity $\vec{v}$ is parallel to the magnetic field, the force is zero and the particle continues in a straight line. If the initial velocity is only in the plane perpendicular to the magnetic field, the particle will travel in a circle. If the initial velocity has a component in the plane perpendicular to the magnetic field and a component parallel to the field, the particle will follow a helix and spiral around the direction of the magnetic field.</p>

<form name="ball">
<p>$m=$ <input type="text" size="20" value="9.1093897E-31" name="mball"> kg $q=$ <input type="text" size="20" value="-1.60217733E-19" name="qball"> C</p>
<p>$\vec{B}=$ <input type='text' size='10' name='Bx' value="0"> $\hat{x} + $ <input type='text' size='10' name='By' value="0"> $\hat{y} + $ <input type='text' size='10' name='Bz' value="5"> $\hat{z}$ [T]</p>

<p style="text-align:center">
<input type="button" value="Load these equations into the form below" onclick="Bfield();">
</p>
</form>

<h3 style="text-align:center">Magnetic deflection of an electron beam</h3>

<p>A magnetic field placed perpendicular to the velocity of electrons in an electron beam will bend the electrons in a direction that is perpendicular to both the velocity and the magnetic field. In the diagram below, electrons are accelerated through a voltage $V_x= 5000$ V towards a positively charged plate. Some of the electrons pass through a small hole in the plate and form and electron beam that travels to a circular region where a magnetic field $\vec{B} = 0.0025\,\hat{z}$ T is generated by a solenoid. The Lorentz force deflects the electrons and they follow a circular path in the magnetic field. The beam then travels further and strikes a screen that is 15 cm to the right of the center of the solenoid.</p>

<p style="text-align:center"><img src="../magnetism/magnetic_deflection.png"></p>

<p>A simulation of magnetic deflection can be found at <a href="http://lampz.tugraz.at/~hadley/physikm/script/magnetism/magnetic_deflection_F.en.php">Magnetic deflection of an electron beam</a>.</p>

<h3 style="text-align:center">J. J. Thomson's experiment to determine the charge-to-mass ratio of electrons</h3>

<p><a href="https://en.wikipedia.org/wiki/J._J._Thomson">J. J. Thomson</a> performed experiments to show that atoms consisted of sub atomic particles that had positive and negative charges. He determined that the negatively charged particles (electrons) were much lighter than the positively charged particles.  The electron charge-to-mass ratio was measured by accelerating the electrons through a voltage $V_x$ towards a positively charged plate. Some of the electrons pass through a small hole in the plate and form and electron beam that travels to a region where an electric field and a magnetic field was present. The electrons were accelerated in the $x$-direction, the electric field was in the $y$-direction and the magnetic field was in the $z$-direction.</p>

<p style="text-align:center"><img src="../magnetism/thomson.png"></p>

<p>The force on an electron due to the electric field is $\vec{F}_{\text{elec}}=-|e|E_y\hat{y}$ and the force on an electron due to the magnetic field is $\vec{F}_{\text{magn}}=|e|v_xB\hat{y}$. The velocity of the electrons can be determined because as the electrons are accelerated in the $x$-direction, the potential energy $|e|V_x$ is converted into kinetic energy, $\frac{1}{2}mv_x^2=|e|V_x$. Therefore $v_x=\sqrt{2|e|V_x/m}$. If the experiment is adjusted so that the electric and magnetic forces cancel each other out $|e|v_xB\hat{y}-|e|E_y\hat{y}=0$, then the charge-to-mass ratio is given by,  </p>

<p style="text-align:center">$\large \frac{e}{m}=\frac{V_y^2}{2V_x\mu_0^2n^2I^2d^2}.$</p>

<p>The calculation below is for $V_x = 5000$ V, $V_y = 60$ V, $\vec{B}= 0.00025\,\hat{z}$ T. Initially, the electrons have zero velocity. The metal plates are 10 cm long with a spacing $d=$ 2 cm. More details can be found at <a href="http://lampz.tugraz.at/~hadley/physikm/script/magnetism/thomson_F.en.php">J. J. Thomson's experiment to determine the charge-to-mass ratio of electrons</a>.</p>

<h3 style="text-align:center">An electron moves in the magnetic field caused by a current carrying wire. </h3>

<p>A magnetic field is generated by a current flowing through a long straight wire in the +$z$ direction,</p>

<p style="text-align:center">$\vec{B}(\vec{r})=\frac{\mu_0I}{2\pi(x^2+y^2)}\,\hat{z}\times (x\hat{x} + y\hat{y})\,\,\text{[T]}$</p>

<p>where $I$ is the current. For details see: <a href="http://lampz.tugraz.at/~hadley/physikm/script/magnetism/line.en.php">The magnetic field around an infinitely-long straight wire</a>.</p>

<p>The force on a particle with charge $q$ in this magnetic field is,</p>

<p style="text-align:center">$\vec{F} =q\vec{v}\times\vec{B}.$</p>

<p>The values for an are initally loaded into the form. The electron moves along a complicated trajectory around the wire. Notice that the work performed (the last column in the table below) remains zero (within numrical error) because the force is always perpendicular to the velocity.  </p>

<form name="elec">
<p>$m=$ <input type="text" size="20" value="9.1093897E-31" name="mball"> kg $q=$ <input type="text" size="20" value="-1.60217733E-19" name="q"> C $I=$ <input type="text" size="20" value="1" name="I"> A</p>

<p style="text-align:center">
<input type="button" value="Load these equations into the form below" onclick="Bwire();">
</p>
</form>

<h3 style="text-align:center">A charged particle in electric and magnetic fields</h3>

<p>When a charged particle with a charge $q$ and a mass $m$ moves in a electric field $\vec{E}$ and a  magnetic field $\vec{B}$, the force on the particle is,</p>

<p style="text-align:center">$ \vec{F} = q(\vec{E} + \vec{v}\times \vec{B}).$</p>

<p>Written out in terms of its three components, the Lorentz force is,</p>

<p style="text-align:center">
$F_x = q(E_x+v_yB_z-v_zB_y)$, <br />$F_y = q(E_y+v_zB_x-v_xB_z)$, <br /> $F_z = q(E_z+v_xB_y-v_yB_x).$
</p>

<p>The components of the electric and magnetic fields can be specified in the text boxes below these can then be loaded into the general 3-D motion differential equation solver to calculate the trajectory of the charged particle.</p>

<form name="elec">
<table>
<tr><td>  $m=$ <input type="text" size="10" value="9.11E-31" name="m"> kg</td><td> $q=$ <input type="text" size="10" value="-1.6022E-19" name="q"> C</td></tr>
<tr><td>  $E_x(x,y,z,t)=$ <input style="border: solid 1px blue" type="text" size="20" value="1E5*sin(1E8*t)" name="Ex"> V/m</td><td> $E_y(x,y,z,t)=$ <input style="border: solid 1px blue" type="text" size="20" value="0" name="Ey"> V/m</td><td> $E_z(x,y,z,t)=$ <input style="border: solid 1px blue" type="text" size="20" value="0" name="Ez"> V/m</td></tr>
<tr><td>  $B_x(x,y,z,t)=$ <input style="border: solid 1px blue" type="text" size="20" value="0" name="Bx"> T</td><td> $B_y(x,y,z,t)=$ <input style="border: solid 1px blue" type="text" size="20" value="0" name="By"> T</td><td> $B_z(x,y,z,t)=$ <input style="border: solid 1px blue" type="text" size="20" value="exp(-x*x-y*y-z*z)" name="Bz"> T</td></tr>
</table>

<p style="text-align:center"><input type="button" value="Load these equations into the form below" onclick='EBfield();'>
</p>
</form>

<p>When a constant uniform electric field is perpendicular to a constant uniform magnetic field, the average velocity of the charged particle will be in the direction perpendicular to both the electric field and the magnetic field. Note that electrons move very fast so a short time step $\Delta t$ must be selected.</p>

<h3 style="text-align:center">Mass-nonlinear spring system</h3>

<p>For a linear spring, the force is proportional to the extension of the spring, $F=-kx$. For a hard spring, the force grows faster than linearly as the spring is extended. An example of a hard spring force is $F=-kx|x|$. For a soft spring, the force grows slower than linearly as the spring is extended. An example of a soft spring force is $F=-k\sin(x)$. If friction is neglected, the mass oscillates around the equilibrium position of the spring. For nonlinear springs, the oscillation frequency depends on the amplitude of the oscillations. For hard springs, the frequency increases with amplitude and for soft springs, the frequency decreases with amplitude. </p>

<form name="ball">

<p style="text-align:center"><input type="radio" id="hard" name="spring" value="hard" onclick="nspring();" checked>
<label for="hard">Hard spring $F=-kx|x|$</label>  <input type="radio" id="lin" name="spring" value="lin" onclick="nspring();">
<label for="lin">Linear spring $F=-kx$</label>  <input type="radio" id="soft" name="spring" value="soft" onclick="nspring();">
<label for="soft">Soft spring $F=-k\sin(x)$</label></p>

<p style="text-align:center">$m=$ <input type="text" size="20" value="0.4" name="mball"> kg $k=$ <input type="text" size="20" value="20" name="kball"> N/m</p>
</p>

<p style="text-align:center">
<input type="button" value="Load these equations into the form below" onclick="nspring();">
</p>
</form>

<h3 style="text-align:center">Resonance: Damped driven mass-spring system</h3>

<p>A mass $m$ is attached to a linear spring with a spring constant $k$. A drag force acts on the mass that is in the opposite direction as the velocity $F_{\text{drag}}=-bv_x$ where $b$ is the drag force constant. The system is driven by a periodic drive force $F_0\cos(\Omega t)$. The differential equation that describes the motion is,</p>

<p style="text-align:center">$\large m\frac{d^2x}{dt^2}+b\frac{dx}{dt}+kx=F_0\cos(\Omega t).$</p>

<form name="ball">

<table align="center">

<tr><td><p style="text-align:right; margin-right:0cm;">$m=$ <span id="mtext">1</span> [kg]</td><td><input type="button" value="-" onclick="ball.slide_m.value=ball.slide_m.value-0.1; resonance();"><input type='range' min='0.1' max='2' step='0.1' name='slide_m' value='1' onchange='resonance();' oninput='resonance();'/><input type="button" value="+" onclick="ball.slide_m.value=eval(ball.slide_m.value) + 0.1; resonance();"></td></tr>

<tr><td><p style="text-align:right; margin-right:0cm;">$b=$ <span id="btext">0.2</span> [kg/s]</td><td><input type="button" value="-" onclick="ball.slide_b.value=ball.slide_b.value-0.01; resonance();"><input type='range' min='0' max='1' step='0.01' name='slide_b' value='0.2' onchange='resonance();' oninput='resonance();'/><input type="button" value="+" onclick="ball.slide_b.value=eval(ball.slide_b.value) + 0.01; resonance();"></td></tr>

<tr><td><p style="text-align:right; margin-right:0cm;">$k=$ <span id="ktext">0.9</span> [N/m]</td><td><input type="button" value="-" onclick="ball.slide_k.value=ball.slide_k.value-0.01; resonance();"><input type='range' min='0' max='1' step='0.01' name='slide_k' value='0.9' onchange='resonance();' oninput='resonance();'/><input type="button" value="+" onclick="ball.slide_k.value=eval(ball.slide_k.value) + 0.01; resonance();"></td></tr>

<tr><td><p style="text-align:right; margin-right:0cm;">$F_0=$ <span id="Atext">0.1</span> [N]</td><td><input type="button" value="-" onclick="ball.slide_A.value=ball.slide_A.value-0.01; resonance();"><input type='range' min='0' max='1' step='0.01' name='slide_A' value='0.1' onchange='resonance();' oninput='resonance();'/><input type="button" value="+" onclick="ball.slide_A.value=eval(ball.slide_A.value) + 0.01; resonance();"></td></tr>

<tr><td><p style="text-align:right; margin-right:0cm;">$\Omega=$ <span id="omgtext">1</span> [rad/s]</td><td><input type="button" value="-" onclick="ball.slide_omg.value=ball.slide_omg.value-0.1; resonance();"><input type='range' min='0' max='10' step='0.1' name='slide_omg' value='1' onchange='resonance();' oninput='resonance();'/><input type="button" value="+" onclick="ball.slide_omg.value=eval(ball.slide_omg.value) + 0.1; resonance();"></td></tr>

</table>

</form>

<p><span id="termv">The resonance frequency is $\omega=\sqrt{k/m-b^2/4m^2}=$ 0.943 rad/s. The maximum response will be when the drive frequency equals the resonance frequency. The amplitude of the steady-state oscillations are $F_0/\sqrt{(k-m\Omega^2)^2+\Omega^2b^2}=$ 0.447 m.</span></p>

<p style="text-align:center"><input type="button" value="Pure resonance" onclick="ball.slide_m.value =1; ball.slide_b.value =0; ball.slide_k.value = 1; ball.slide_A.value = 0.1; ball.slide_omg.value = 1; ode.Nstep.value = 3000; resonance();" > <input type="button" value="Underdamped" onclick="ball.slide_m.value =1; ball.slide_b.value =0.1; ball.slide_k.value = 1; ball.slide_A.value = 0.1; ball.slide_omg.value = 1; ode.Nstep.value = 3000; resonance();" > <input type="button" value="Critically damped" onclick="ball.slide_m.value =0.5; ball.slide_b.value =1; ball.slide_k.value = 0.5; ball.slide_A.value = 0.1; ball.slide_omg.value = 1; ode.Nstep.value = 3000; resonance();" > <input type="button" value="Overdamped" onclick="ball.slide_m.value =0.5; ball.slide_b.value =1; ball.slide_k.value = 0.2; ball.slide_A.value = 0.1; ball.slide_omg.value = 1; ode.Nstep.value = 3000; resonance();" ></p>

<p>For more details see: <a href="http://lampz.tugraz.at/~hadley/physikm/script/oscillations/resonance.en.php">Resonance of a damped driven harmonic oscillator</a>.</p>

<h3 style="text-align:center">Swinging spring</h3>

<p>The swinging spring or the <a href="https://en.wikipedia.org/wiki/Elastic_pendulum">elastic pendulum</a> is a mass connected to one end of an elastic rod that stretch but cannot bend. The other end of the rod is attached at one fixed point so that the mass can swing freely. This system can execute complex motion and for some initial conditions it executes chaotic motion. The equations are:</p>

<p style="text-align:center">$\large m\frac{d^2x}{dt^2}=-m\omega^2\left(\frac{|\vec{r}| -l}{|\vec{r}|}\right)x$</p>
<p style="text-align:center">$\large m\frac{d^2y}{dt^2}=-m\omega^2\left(\frac{|\vec{r}| -l}{|\vec{r}|}\right)y$</p>
<p style="text-align:center">$\large m\frac{d^2z}{dt^2}=-m\omega^2\left(\frac{|\vec{r}| -l}{|\vec{r}|}\right)z-mg$</p>

<p>where $m$ is the mass of the rod,  $l$ is the unstretched length of the rod,  $g = 9.807$ is the acceleration of gravity at the earth's surface, $\omega^2=k/m$, and $|\vec{r}| =\sqrt{x^2+y^2+z^2}$. </p>

<form name="swgspr">
<table>
<tr><td>  $m=$ <input type="text" size="10" value="1" name="m"> kg</td><td> $k=$ <input type="text" size="10" value="5" name="k"> N/m</td><td>  $l=$ <input \ type="text" size="10" value="1" name="ell"> m</td><td></td></tr>
</table>

<p style="text-align:center"><input type="button" value="Load these equations into the form below" onclick="swingspring();">
</p>
</form>