Mathcad Exercises

 We are moving: This site will move to around September 1, 2012.

Send corrections to scott.robertson [at]

For a summary of the conversion to Prime 2.0 and related issues, scroll down to the area below all of the exercises.

Table of Contents

Preliminaries:  Mathcad Primer

Introduction to Mathcad, written by Prof. Dubson

Short introduction to Mathcad.  This document is NOT a "live" Mathcad document, but rather a printable PDF.  It is still handy, though.

The intro was not written for use with Mathcad Prime. 

Ch.0  Mathcad Essentials

Variables types and simple operations

Examples are given of simple mathematical operations and a simple graph.





This exercise demonstrates how to enter matrices and the matrix operations from the matrix toolbar.

Also illustrated are the rows, cols, stack, augment, reverse and submatrix functions.




Vectors, lookups, and interpolation

This exercise demonstrates the standard vector operations and operations that might be used on data stored as vectors, such as reverse, value lookup, index lookup (match), max, min, linear interpolation, spline interpolation, and extrapolation. 




Boolean functions and If(,,)

This exercise demonstrates the if(,,) function and some of the Boolean functions such as >, <, or, and, not, and xor. Typing hints are given for if(,,).




Plots of vectors, matrices and bitmaps

Examples are given of contour and mesh plots for data that are a function of two variables. Also illustrated is Picture from the matrix toolbar.




Creating and editing programs and loops

This exercise illustrates how to type in programs containing loops. The confusing task of opening new lines within, above, or below loops is illustrated.





This short exercise demonstrates using units and making units conversions.



Ch.1  Electrostatics

Plotting the field of an electrostatic dipole


This exercise begins by using the vector field plot to display the field of a single charge and of a pair of charges.

For the dipole, the lines of force are followed by integrating a short differential distance ds that is parallel to the E field.

Euler's method is used for the integration. Field line plots are also made in the Magnetostatics Exercises below.





Laplace Solver


This exercise shows how successive approximations are used to obtain solutions to Laplace’s equation on a 2-d grid.

Physics/math content: teaches that a characteristic of Laplace’s equation is that the potential at a grid point is the average of the values at the neighboring grid points.

Mathcad content: program loop, if function, relax function, matrices, matrix transpose.

Mathcad graphics: Contour plots





Laplace's Equation with a Dielectric


This exercise solves Laplace's equation first in one dimension for a dielectric slab between biased plates.

Then Laplace's equation is solved on a two dimensional grid for a dielectric rod between biased plates.

The content is similar to that for the Laplace Solver exercise.





Poisson’s Equation and Debye Shielding

This exercise shows how successive approximations are used to obtain solutions to Poisson’s equation in one-dimension.

It uses electron and ion densities determined by the Boltzmann factor.

Physics content: shows that the potential around an electrode placed in plasma decays exponentially with distance. Demonstrates the utility of dimensionless units.

Math content: shows a numerical instability and how to stop it.

Mathcad content: Program loop.

Mathcad graphics: Simple x,y plots.





Poisson's equation in cylindrical geometry


The finite-difference form of Poisson’s equation in cylindrical geometry is solved by the relaxation method, first assuming no z dependence and then allowing a z dependence.

A cylinder of charge is assumed to be within a cylindrical container at zero potential. This extends the previous exercise to cylindrical coordinates.





The ionosphere: Separation of ions by gravity


The ionosphere contains ions of both oxygen, O+, and hydrogen, H+. The proportion of H+ increases with altitude because the H+ is lighter.

The electrons are the lightest species, but are constrained by quasineutrality to be at the same altitude as the ions.

In this exercise we use the relaxation methods for Poisson's equation to apply the quasineutrality constraint and find the altitude dependence of the densities of O+, H+, and electrons.




Ch.2  Magnetostatics

Magnetic Field of Parallel Wires

This exercise plots the field of a single wire, two wires with parallel currents, and two wires with antiparallel currents.

A magnetic “x point” is illustrated. The vector fields are added.

Physics content: The azimuthal field of a wire is projected onto the x and y directions.

Mathcad graphics: Plots a vector field in two dimensions.





Plotting Field Lines


This exercise shows how to integrate with Runge Kutta a unit vector to obtain a curve in 3-d space. The components of Bq are projected onto the x and y directions.

Mathcad content: Demonstrates the Runge-Kutta solution of a differential equation using rkfixed.

Mathcad graphics: illustrates a 3-d scatter plot that can be tilted by the user.





Field of a Loop


This exercise plots the field of a loop of wire. The expression for B involves the complete elliptic integrals E(k) and K(k). These integrals are defined at the outset.

The unit vector is found in the r,z plane and integrated. The result is B looping around the minor axis for points close to this axis.

In a second plot, a constant field is added which makes an x point.

Mathcad content: integration of a unit vector as in the field line plotting exercise. Evaluation of elliptic integrals. TOL is used to decrease calculation time.

Mathcad graphics: 2-d plots of lines.





B Field of the tokamak


This exercise plots the magnetic surface mapped by a magnetic field line of a tokamak.

It follows a magnetic field line created by a current loop perpendicular to the z axis and a wire along the z axis.

The concepts of magnetic surface and aspect ratio are discussed. The content is similar to that in Field of a Loop.




Ch.3  Lorentz Equations of Motion

Gyrating motion: a video

This exercise integrates the Lorentz equations of motion (v x B only) in 2-d with no E field.

There are two coupled equations for vx and vy which are put into a video. A simple plot shows that the two velocities are out of phase.

The effect of time step on error is discussed briefly.

Mathcad content: illustrates Runge Kutta (rkfixed).

Mathcad graphics: demonstrates making a video using the Animate command.

Download video: Gyration4.avi.



Not supported in Prime 2.0


Deflection of a particle beam by a bending magnet


This exercise follows a charged particle beam through a bending magnet. The beam is also focused.

The Runge-Kutta integrator is used to integrate the equations of motion in two dimensions.

Mathcad content: rkfixed.





E x B drift in uniform fields


This exercise integrates the Lorentz equations of motion in 3-d, though two would do for the E x B drift.

There are six equations to be integrated and six phase space variables. The E x B drift is plotted and simple acceleration along B.

Mathcad content: rkfixed. ORIGIN = 1 is used to make subscripts 1,2,3 correspond to x,y,z.





Drifts in the field of a wire


This demonstrates the gradient drift and the curvature drift in the azimuthal field around a straight wire.

The initial conditions are adjusted so the two drifts appear together and appear separately.

This is more sophisticated than the previous exercise because B is now a function of position and must be calculated.

Mathcad content: rkfixed, vector cross product, ORIGIN.

Mathcad graphics: a 3-d scatter plot with gradient shading and with two arguments. ORIGIN.





Ion motion in the Earth’s magnetic field: Van Allen’s belts

The earth’s field is represented as a dipole and a dipole expansion is used for Bx, By and Bz. There are three vector components of B, which vary spatially.

First, a plot of some field lines from pole to pole. Then an ion is followed. Real units are used so distances are in meters and time in seconds.

Physics content: Demonstrates magnetic mirroring.  

Mathcad content: rkfixed, vector cross product, ORIGIN, stack for matrices.





Particle Orbits in the Multidipole Device


This exercise finds the trajectory of a 1 eV electron in a common type of laboratory plasma device in which a cold plasma is confined by a multidipolar field.

The magnetic field components are a function of only x and y, thus the field values are calculated in advance and placed in a matrix.

During the trajectory calculation, the field values at the electron location are found by a two-d interpolation.

A wire array is used to generate a magnetic field similar to that of the permanent magnets. The magnetic vector potential is calculated and plotted.





Banana orbits in tokamaks

This exercise uses the drift approximation to find the path followed by an energetic ion in a tokamak.

The equation of motion parallel to B for the bounce motion caused by the mirror force is solved simultaneously with the equation for the vertical drift.




Ch.4  Fluid-Poisson System of Equations

Child’s law

The flow of electrons from a cathode to an anode is followed with the potential found from the fluid continuity equation and Poisson’s equation.

The analytic solution is compared with the Runge Kutta solution.

Mathcad content: rkfixed. A method is given for solving a differential equation for a range of values of a coefficient in the equation.

The “trick” is to put the coefficient in the list of boundary conditions.





Bohm’s sheath

The fluid continuity equation and Poisson’s equation are solved for the sheath that occurs at the boundary of a plasma.

This is similar to Child’s law but has both electrons and ions.

Mathcad content: rkfixed. “Shooting” is demonstrated; root is used to find a starting condition by trial and error. 





Sheath and presheath

The potential profile from the midplane to the wall is found using Poisson’s equation.

The electron density is modeled by the Boltzmann relation and the ion density is found by solving the continuity and momentum equations for the ions.

A homogeneous source of ionization is assumed.




Ch.5  Waves

Electrostatic waves in a cold magnetized plasma

The dispersion relation for a magnetized cold plasma with electrons and ions is solved.

For a given k and angle of propagation, a root finder is used to find the frequencies of waves.

The upper and lower hybrid frequencies are found for perpendicular propagation and the plasma frequency is found for parallel propagation.

Physics content: frequencies of waves in magnetized plasma

Mathcad content: root finding on an interval or with an initial guess, using root.





Extraordinary wave dispersion relation

Similar to the above exercise except this wave is electromagnetic. The vectors E and k are perpendicular to B.

There are two roots to the dispersion relation that are found using the polyroots function.

Mathcad content: polyroots





Ray tracing in an inhomogeneous medium

The ray of a radio wave is followed as it approaches the ionosphere. The ray bends gradually and descends back to the earth’s surface.

The differential equations that advance the position x and the wavevector k are presented and solved by Runge Kutta in three dimensions.

The equations are based on the WKB approximation. The dispersion relation w(k,x) is a function of altitude because of the increasing electron density.  5_Waves_Ray_tracing.xmcd



Ch.6  Waves with Complex w

Two stream instability

This looks at the fluid dispersion relation for a beam of cold electrons flowing through a background of ions.

Mathcad content: uses the root finder root to search for complex roots of an expression given a guess. Uses polyroots for find the roots of a polynomial expression.

Graphics content: simple x, y plots.





Landau dispersion relation and Landau damping

This exercise evaluates the dispersion relation for electrostatic plasma waves using the Vlasov equation.

The real part of the dispersion relation is expressed in terms of the W function.

Both the principal value integral is found and the value from expansion of the denominator as a power series.

The imaginary part is plotted showing that waves with k < 0.2 are weakly damped.

Math content: principal value integral.

Mathcad content: Evaluation of integrals with integrable singularities.




Ch.7  Random Walks, Collisions, Resistivity and Diffusion 

Random numbers and random events

Uses the rnd(1) function to generate random numbers on the interval 0,1.

A logarithmic “stretch” is used to generate the time intervals between random events. The decay of a population of 1000 nuclei is simulated.

Lastly, speeds are selected from a Maxwellian distribution by the rejection method that is applicable to any distribution function.

Mathcad content: rnd, hist, if, program loop.

Mathcad graphics: Bar charts, stem charts and histogram function.





Random walks and diffusion

The trajectories of particles are followed with randomly occurring collisions that change the direction of the velocity vector but not the length.

This corresponds to elastic electron-neutral collisions. A collection of particles is started at the origin and the spreading of the locations with time is followed.

This is done both in two and three dimensions, with random walks characterized either by a mean time between collisions or a mean free path between collisions.

The mean square distance is compared with that from the diffusion equation.

Physics content: collisions with constant mean free path or with constant mean free time.

Mathcad content: same as in Random numbers and random events.






An electric field is added to the random walk above. The particles then attain a drift velocity from which the current and resistivity can be determined.

The mean squared velocity increases with time showing resistive heating. It is shown that the results conserve energy.

Physics content: resistivity, mobility, resistive heating.





Diffusion equation

The diffusion equation is solved using the simplest finite-difference method in Numerical Recipes.

A top hat initial distribution is evolved in time showing the rapid rounding of the corners (higher order modes) and slower decay of the lowest order mode.

The decay of the lowest order mode is shown to agree with the solution obtained by separation of variables.

Math content: a finite difference method for the diffusion equation.

Mathcad content: Matrix operations. A matrix defined as a function of variables. Program loop.

Mathcad graphics: 3-d surface plot with time on one axis.





Diffusion in Cylindrical Geometry


This simply repeats the Diffusion Equation exercise in cylindrical geometry.




Ch.8  Kinetic Theory

Liouville’s equation

The flow of electrons is followed from a thermionic cathode to an anode as in Child’s law. The initial distribution function is a half-Maxwellian and the evolution of this distribution is followed showing that an accelerating beam becomes “colder” as a consequence of the preservation of phase space density.





Fokker-Planck Equation

The Fokker-Planck Equation describes the evolution of a distribution function as a consequence of collisions of the particles among themselves.

This exercise applies the Lax-Wendroff method to solving the Fokker-Planck equation.

In the first example, a distribution function with a cutoff tail is followed in time, showing that the tail is filled-in by collisions.

In the second example, a beam of particles is followed as it is slowed by collisions.





Landau Damping from a Numerical Solution to the Vlasov Equation

The Cheng-Knorr algorithm is used to solve the Vlasov-Poisson system of equations to find Landau damping of an electrostatic wave.





Source and Collector Sheaths

This extends the treatment of the sheath problem to a plasma of two species. The Q machine has a hot plate that emits both electrons and ions.

This exercise finds the potential profile in the device (1) with emitters at both ands and (2) with an emitter at one end and a collector at the other end.

The collector is assumed floating and the potential on the collector is updated continually to be consistent with the plasma properties.

Kinetic theory is used to find the left-going and right-going electron and ions distribution functions at each grid point.





Sheath Above an Electron Emitter

A sheath composed only of electrons, a "virtual cathode", forms above a surface that emits electrons, such as a thermionic or photoelectric emitter.

This exercise uses Poisson's equation to find the potential profile between parallel planes, one of which is emitting a one-sided Maxwellian distribution of electrons. Kinetic theory is used to find the electron distribution at each point between the plates.




Ch.9  Data Analysis

Langmuir Probe Data Analysis

This exercise analyzes the data from a cylindrical Langmuir probe and finds the electron temperature and electron density.

A second more sophisticated analysis resolves the electron distribution into the sum of two populations with different temperatures.





Druyvesteyn's Probe Analysis

In 1930 Druyvesteyn showed that the electron speed distribution is proportional to the second derivative of the probe data.

This lengthy exercise imports probe data, removes the ion contribution to the probe current, and finds the electron speed distribution.

For comparison, it also finds the two Maxwellian distributions that best describe the data.

A new method is used to find the distribution function projected on the the plane perpendicular to the probe axis. This method is less affected by noise.




Ch.10  Miscellaneous Exercises

Harmonic oscillator

This exercise introduces Runge-Kutta to solve equations of motion in one dimension. The second order differential equation is expressed as two first order equations.

Mathcad content: illustrates Runge Kutta (rkfixed).





Harmonic Oscillator with Energy Conservation

This exercise is a second way of solving  the harmonic oscillator problem that uses conservation of energy.

The velocity is found as a function of position using energy conservation, and the velocity is integrated to find the position as a function of time. Energy conservation determines the square of the velocity, and whether the positive or negative value of the square root is used is determined using the sign of the velocity obtained from a linear extrapolation.

Mathcad content: illustrates Runge Kutta (rkfixed).





Orbital motion

An electron orbits an ion classically. This could also be celestial mechanics. Turning on a weak perturbation destroys the spherical symmetry and causes the orbit to tilt.

This demonstrates the electric part of the Lorentz force in three dimensions. Introduces using Runge Kutta for three-dimensional motion. A six dimensional phase space is introduced.

Mathcad content: illustrates Runge Kutta (rkfixed), and the stack and sub matrix commands.

Graphics content: 2-d and 3-d point plots.





Paul trap, ponderomotive force, and Mathieu’s equation

The motion of an ion in a Paul trap is followed in the r,z plane. Containment is from the ponderomotive force exerted by AC potentials on external electrodes.

The problem is in (r,z) coordinates and the equations of motion along r and z are solved by Runge Kutta.

The relation to Mathieu’s equation is discussed and the student is encouraged to examine stability.

Physics content: ponderomotive force.





Standard Map

This exercise examines the Standard Map of nonlinear dynamics. A program loop is used to find the mappings for a wide variety of seed values.





Fluid flow from a potential function (new in 2008)

The relaxation method is used to find a potential function that describes the incompressible flow of fluid in a nozzle.

A second solution is found with an object placed in the flow.




Fluid flow from a stream function (new in 2008)

The relaxation method is used to find a stream function that describes the incompressible flow of fluid in a nozzle.

A second solution is found that is circular flow about an object placed in the flow.




Equilibria of a Cylindrical Plasma (new in 2008)


The MHD equilibrium equation J x B - grad P = 0 is integrated to find the current and magnetic field profiles in cylindrical geometry:

First with pressure and no axial current, second with axial current and no pressure, and third with both axial current and pressure.




Axisymmetric equilibria from the Grad-Shafranov equation (new in 2008)

The Grad-Shafranov equation is derived and used to find equilibria satisfying J x B - grad P = 0 for a tokamak.

A relaxation method is used that is similar to the method used for Poisson's equation.

The program finds the current densities Jr, Jz, and Jf as well as the fields Br and Bz.

The inputs are the rectangular boundaries, the pressure on the magnetic axis, the safety factor q on the magnetic axis, and the applied toroidal field Bf




Ch.11  Finite Difference Methods

Numerical Integration of a Function

This exercise introduces the concept of a grid of points and finds the area under a curve by Simpson’s rule.





Numerical Integration of Derivatives

This exercise compares integration by Euler’s method, by the midpoint method, and by the fourth order Runge-Kutta method.





Solving a Second Order Differential Equation

This exercise demonstrates changing a second-order differential equation into two first order equations that can be solved by the Runge-Kutta method.





Derivation of Finite Difference Equations

This exercise uses the Taylor expansion to derive finite-difference methods for solving differential equations. The Runge-Kutta method is derived for the special case of a function depending only upon the independent variable.




Ch.12  Circuit Analysis

Circuit Analysis with Complex Variables

In this exercise, we use complex variables for circuit impedances and find the response function for RC low pass, RC high pass, and RLC bandpass filters.

Both the amplitude and the phase of the signals are plotted as a function of frequency.





Response of an RC Filter to a Square Wave

In the lab, square waves or timing pulses are often passed through circuits that filter the wave.

In this exercise, we find what a square pulse will look like after passing through a low-pass, a high-pass, and an RLC bandpass filter.

The square wave is represented by sum of Fourier components.





Reducing Noise by Signal Averaging

A signal that is buried in random noise can often be recovered by signal averaging methods.

In this exercise we generate a vector of digitized voltages similar to the vector that would be generated by a digital oscilloscope that is triggered by a synchronizing pulse.

We then show that averaging signals improves the signal-to-noise ratio.




Ch.13  Particle-in-cell (PIC) methods for plasmas (new in 2011)

The particle-in-cell method was developed to study electrostatic waves, beam-plasma interactions, sheaths, and other phenomena in plasma that require simultaneous solutions to Poisson's equation and the equations of motion for electrons and (sometimes) ions. In the exercises below PIC codes are developed for problems with increasing levels of complexity. 


Basic Particle-in-Cell Code for Cold Plasma Oscillations (new in 2011)

This basic code solves Poisson's equation simultaneously with the equations of motion for electrons. An electrostatic wave is made by modulating the electron velocities.

Simplifications include: (1) the ions are assumed to be a uniform neutralizing background and (2) the electrons have no thermal velocity so they do not leave the simulation domain.




Particle-in-Cell Code for a Cathode Sheath (new in 2011)

This code finds the electrostatic sheath between an emissive cathode and an anode in a vacuum tube. There are no ions.

The code keeps track of the number of electrons in the simulation domain which varies as electrons are emitted from the cathode, return to the cathode, or are lost at the anode.

The anode may have a bias potential.





Particle-in-Cell Code for Landau Damping (new in 2011)

This code differs from the one for cold plasma oscillations in that the electrons are given initial velocities chosen from a distribution function.

Periodic boundary conditions are assumed and electrons leaving the simulation domain at one boundary are returned at the other boundary.

An electrostatic wave is launched by modulating the electron density and this wave is observed to damp at approximately the predicted Landau Damping rate.

Plasma wave echoes are observed and the damping is absent for a "top hat" distribution.




Particle-in-Cell Code for Two-Stream Instability (new in 2011)

This code is similar to the Landau Damping code except that the electron distribution function is two oppositely directed streams each described by a Maxwellian distribution with a drift velocity.




Particle-in-Cell Code for Beam Plasma Instability (new in 2011) 

This code is similar to the Landau damping code except that there is an electron beam in addition to the Maxwellian plasma which creates a “bump-on-tail” instability.




Mathcad Prime 2.0

·         We have added all the programs we converted to Mathcad's new platform - Prime 2.0.

·         The good: no calculation errors were encountered when checking the converted worksheets.

·         The bad: you can’t do the conversion with your existing version of Mathcad – even Mathcad 15. PTC (the software company) makes you download a 30-day trial version of Mathcad 15 in order to use the converter in Prime. This trial version of Mathcad 15 is version M010. The converter in Prime works ONLY with Mathcad 15 version M010 (the aforementioned trial version).

·         The ugly: installing the Mathcad 15 trial version removes any prior version of Mathcad you might have had on your computer. Once the 30-day trial period expires, you lose the ability to use the old Mathcad and can only use Prime from that point forward. Therefore, you must convert any old Mathcad doc(s) you want within that timeframe. Although nothing happens to your old files if you don’t convert them, once the 30-day trial elapses you will have to again install the Mathcad 15 trial version to perform any conversions. In addition, functionality decreased significantly in Prime and many features are no longer supported. This is explained more thoroughly below.


·         Prime 2.0 Documentation: here are useful pdf files







As mentioned in the above pdf files, 'annotations' will appear in a converted worksheet to let you know if there are possible issues related to the conversion. The vast majority of these annotations just point out places where a formatting change was made and can be ignored. However, Prime only annotates a relatively short list of possible issues - many other issues were discovered by going through the exercises. The exercises for chapters 1 through 13 were scrubbed before being posted. Here is what we have so far:

1.      Data and calculations were checked first - formatting might be slightly different but that is all.

2.      Next, we attempted to restore the original layout of the worksheets with regard to content which might have migrated to an adjacent or subsequent page. Prime does not have a print preview option, so be careful before printing anything.

3.      The page numbers and date in the header don't convert. The header box is positioned so near the top of the page that most printers will not print it. The header should be 0.5 in from the top of the page and is not.

4.      Text boxes have lost superscripts and subscripts. So if we said 3 x 1010 in a box, that is changed to 3 x 1010. And if we referred to x2 it became x2. And if we talked about a subscripted variable xj, it became xj.

5.      We put some equations from the Microsoft Equation editor in our Mathcad 15. It is not longer possible to edit these or make new ones.

6.      Our videos in Mathcad 15 don't work in Prime.

7.      Our vector plots in Mathcad 15 don't work in Prime. The vector plot option is gone.

8.      Plot labels are not supported. Axis labels are still available but if you want a plot label then you will need to make a new text box that serves as your plot label.

9.      Results can't be displayed as tables anymore - they are converted to a matrix.

10.  Only one contour plot can be plotted at a time.

11.  Point plots with more than 4,998 points are automatically made into line plots.

12.  Plot tick intervals are usually different after conversion, so you have to change them back.

13.  To be continued...

Visual to help with annotations you might come across:

Annotations in Prime.pdf