# Streamlines

Jump to: navigation, search

# 2. Overview

In scientific visualization, a scalar is a single value that is associated with each point, cell, or cell element. The most common type of scalar plot that you may be familiar with is temperature versus position is a region on the surface of Earth . Scalars are typically visualized using contour lines and images with colorbars, both of which are shown in .

Streamlines are commonly used to visualize vectors measurements as an alternative to, or in addition to, drawing many arrows. In the following image, velocities are measured at 100 points (indicated by dots) on a regular grid. The most straightforward way of visualizing the data is to draw arrows with length proportional to the magnitude of the vector. (The program  was used to create the image.) A streamline is the path a particle would follow if placed in the velocity field (assuming that the particle's presence would not modify the velocity field and that the velocity field does not change as the particle moved). The path of a particle can usually be approximately found by looking at the vectors. However, to get an exact answer, calculations are required.

Let px(1) and py(1) be the initial position of the particle. Both px and py are arrays that will have the position of the particle at various instants of time. The subsequent position of the particle will be determined by taking a series of small steps; for a massless particle, the length of the steps in the x- and y-directions is proportional to the velocities in the x- and y-direction, with a proportionality symbol of dt. This means that the particle always has a velocity that is equal to the measured velocity at a given location.

If the particle is released at px(1) = 1 and py(1) = 1, it is expected to initially follow the path of the vector drawn at that position (vx(1)=1.5 vy(1)=2.0). The procedure for determining the location after one step in this direction is

px(2) = px(1) + dt*vx(1)
py(2) = py(1) + dt*vy(1)


where dt is a small number with units of time.

(Note that there exist better algorithms for taking steps with higher accuracy, but this is the easiest algorithm to interpret.) The smaller the value of dt, the smaller the overall step size (and the more accurate the computed path). Assuming dt = 0.1, the position after the first step is

px(2) = 1.0 + 0.1*1.0 = 1.1
py(2) = 1.0 + 0.1*1.5 = 1.15


This position is shown in the following figure, which is a zoomed-in version of the previous figure. To compute the position after the second step, we need to evaluate

px(3) = px(2) + 0.1*vx(2)
py(3) = py(2) + 0.1*vy(2)


Because the new position of the particle px(2),py(2) is not exactly at a point where we have a measurement of velocity (indicated by black dots in the figue), we must approximate the velocity. All four of the nearest measurements to (px(2),py(2)) = (1.1,1.15) have vx = 1.0, so we can assume vx(2) = 1.0. To estimate vy(2), we can use linear interpolation and show that its estimate gives

vy(2) = 1.55


Note that in general, bi-linear interpolation is needed to determine both vx and vy. See #Code Example for an example.

Given the estimates of vx(2) and vy(2), we can compute the predicted third position:

px(3) = px(2) + 0.1*vx(2) = 1.1  + 0.1*1.0  = 1.2
py(3) = py(2) + 0.1*vy(2) = 1.15 + 0.1*1.55 = 1.305


And then these two points can be used to compute the next position

px(4) = px(3) + 0.1*vx(3)
py(4) = py(3) + 0.1*vy(3)


after estimating vx(3) and vy(3) in the same was that was done for vx(2) and vy(2).

The following image shows the path of a massless test particle after taking many steps and using the above procedure to compute the positions. Many such streamlines can be drawn by repeating the above procedure changing only the starting position for the streamline. # 3. Interpolation Review

Note that more details on interpolation are given at Interpolation.

## 3.1. 1-D

You have probably already dealt with linear interpolation: you have to points with an associated scalar value and want to estimate the scalar value somewhere in-between. For example, suppose that you have temperature measurements at two points:

x T
1 0.1
2 0.3


In your math class, you would determine the equation of the line that connected these two points and then use it to determine an estimate of temperature at an intermediate point. This is called linear interpolation or 1-D interpolation. In MATLAB, the calculation of Te, the estimate of T at x = 1.1 is

Te = interp1([1,2],[0.1,0.3],1.1)


The first two arrays are the position and temperature, respectively. MATLAB uses this information to compute an estimate of the temperature at other points, given by the last argument to interp1.

In the example given in #Overview, measurements of velocity were available at four boundary points and we needed to estimate vx and vy at a point in the interior. The velocity measurements on the boundary happened to be such that only 1-D interpolation was needed to estimate vy, because vy only varied in the x-direction.

In this case, the measurement table is

x  vy
1  1.5
2  2.0


and the estimate of vy at x = 1.1 is

vy = interp1([1,2],[1.5,2.0],1.1) % Gives 1.55


## 3.2. 2-D

 Now, suppose that you have a set of measurements of height at four grid points and want to estimate the height somewhere within inside. x y H 1 1 10 1 2 20 2 1 30 2 2 40  The function interp2 can be used to do bi-linear or 2-D linear interpolation. However, MATLAB now requires three inputs: the X positions, the Y positions, and the heights, and it requires them to be given as matrices. In addition, as described below, the matrices must values in a certain order. X = [ 1 , 2  ; 1 , 2] Y = [ 1 , 1  ; 2 , 2] H = [10 , 30  ; 20 , 40] % Estimate height at (x,y) = (1.1,1.1) % X, Y, and H are needed in order to make the estimate. He = interp2(X,Y,H,1.1,1.1) % Answer should be closest to value at (x,y) = (1,1) Important Note

MATLAB requires the arrays X and Y to have a certain form. For example, one may expect that we could re-arrange the order of the values in X and Y provided that we re-ordered the values in H for consistency:

X = [ 1  ,   1  ;  2  ,  2]
Y = [ 1  ,   2  ;  1  ,  2]
H = [10  ,  20  ; 30  , 40]


The above set-up is another way of representing the information in the original table:

x   y   H
1   1   10
1   2   20
2   1   30
2   2   40


However, this

He = interp2(X,Y,H,1.1,1.1)


results in an error. The reason is that interp2 requires the values in the X matrix to be increasing from left to right (when displayed on the screen) and the values in Y to be increasing from top to bottom. In contrast, the example given for interp1 will work even if we reverse the order of both x and T: interp1([2,1],[0.3,0.1],1.1) give the same result as interp1([1,2],[0.1,0.3],1.1).

## 3.3. 2-D and 1-D equivalence

Sometimes a problem that appears to require 2-D interpolation can actually be solved using 1-D interpolation. This can happen if all of the variation of the parameter to be interpolated is in one direction (example given here) or if the point of interest is on the edge of the box with corners given by the x and y values (this second case is demonstrated in the problems section.

Suppose that you have a set of measurements of height and temperature at four grid points as given in the table and want to estimate the height somewhere inside of the four points. A plot of the values in the table is shown to the right.

 x y H 1 1 10 1 2 20 2 1 10 2 2 20  Because there is only variation in H in the y direction, one could use 1-D interpolation to find the value of H - the equivalent problem is to find the height at y=1.1 in y H 1 10 2 20  and note that this height estimate applies to all values of x between 1.0 and 2.0. # 4. Code Example

The following program shows the various steps required to draw a streamline.

First, position matrices X and Y are created. A given element in X and Y represents the x- and y-positions of a measurement. E.g. X(1,1) and Y(1,1) represent a point where a measurement was made.

Next, matrices U and V are created. Each element in U corresponds to an x velocity measurement made at a position given by an element in X and Y. E.g., U(1,1) is the x velocity of a measurement made at position x=X(1,1) and y=Y(1,1). The rest of the steps are given in the code that follows. Each element in V corresponds to a y velocity measurement made at a position given by an element in X and Y.

% Specify matrix of x positions
% interp2 requires that this matrix has values that increase from left to right
X = [1   ,  2   ; 1   , 2]
% Specify matrix of y position
% interp2 requires that this matrix has values that increase from top to bottom
Y = [1   ,  1   ; 2   , 2]
% Specify matrix of velocities in x direction
U = [ 1.0,  1.0 ; 1.0, 1.0]
% Specify matrix of velocities in y direction
% The location of the velocities in this matrix depend on how X and Y are set up.
% For example, the value 3.0 at location V(2,1) is measured at position X(2,1) and Y(2,1).
V = [ 1.0,  2.0 ; 3.0, 4.0 ]

% Code above generates the vectors and their locations. The following code computes streamlines.

px(1) = 1.0;
py(1) = 1.0;
% Here we use the function interp2 to determine an estimate
% of the x velocity at point px(1), py(1).  Because px(1) is exactly on
% a grid point, the estimated velocity should exactly match the velocity given
% in U.
vx(1) = interp2(X,Y,U,px(1),py(1));
% Same as above, except for y velocity
vy(1) = interp2(X,Y,V,px(1),py(1));
% Because px(1) and py(1) correspond to a grid point, interpolation was not needed
% because these positions are at a measurement point.
% We could have just said vx(1) = 1.0 and vy(1) = 1.0, but writing
% with interp2 will allow us to change px(1) and py(1) to be
% a point that is not at a measurement location
% without having to manually change vx(1) and vy(1).

dt = 0.01;
Nsteps = 10;
for t = 2:Nsteps
% Take a step in x.  New position is px(t).  Last position is px(t-1).  Last velocity is vx(t-1).
px(t) = px(t-1) + dt*vx(t-1);

% Take a step in y.  New position is py(t).  Last position is py(t-1).  Last velocity is vy(t-1).
py(t) = py(t-1) + dt*vy(t-1);

% Find velocity at new position.
vx(t) = interp2(X,Y,U,px(t),py(t));
vy(t) = interp2(X,Y,V,px(t),py(t));

% Display some information
fprintf('px(%d) = %.2f  py(%d) = %.2f  vx(%d) = %.2f  vy(%d) = %.2f\n',...
t,px(t),t,py(t),t,vx(t),t,vy(t));

end


# 5. Field Lines

Consider a vector field defined by $\mathbf{F} = F_x(x,y)\hat{x} + F_y(x,y)\hat{y}$

Thus far, we have considered cases where Fx and Fy were external velocities and then noted that at a given instant in time, an object at a given position would have a velocity of

(1) dx / dt = Fx(x,y)

(2) dy / dt = Fy(x,y)

A field line is a line that is drawn from a given point in space such that it is always tangent to the field. To determine the equations that must be solved to find the trajectory, consider a differential distance ds drawn in the direction of $\mathbf{F}$ that is made by stepping a distance dx in the x-direction and dy in the y-direction. This triangle is similar to a triangle formed with Fx in the x-direction, Fy in the y-direction, and F as the hypotenuse so that

dx / ds = Fx(x,y) / | F |

dy / ds = Fy(x,y) / | F |

To determine the x- and y- values that correspond to a line that is always in the direction of $\mathbf{F}$, we can solve the above equations in the same manner as (1) and (2). That is, to compute the positions for (1) and (2), we used the approximations $dx/dt \simeq \Delta x/\Delta t$ and $dy/dt \simeq \Delta y/\Delta t$ and then wrote

x(i+1) = x(i) + dt*Fx(x(i),y(i))
y(i+1) = y(i) + dt*Fy(x(i),y(i))


For computing the field line, we instead use the approximations $dx/ds \simeq \Delta x/\Delta s$ and $dy/ds \simeq \Delta y/\Delta s$ so that

Δx / Δs = Fx(x,y) / | F |

Δy / Δs = Fy(x,y) / | F |

Written out in code form, this is

F = sqrt(Fx(x(i),y(i))^2 + Fy(x(i),y(i))^2)
x(i+1) = x(i) + ds*Fx(x(i),y(i))/F
y(i+1) = y(i) + ds*Fy(x(i),y(i))/F


To solve for x and y that correspond to a line that is always tangent to $\mathbf{F}$, select an inital value for x and y, an appropriate value of ds, and repeat the above three equations for many values of i.

# 6. Problems

## 6.1. 1-D Interpolation

An estimate of temperature at x = 1.1 given the following measurements

x T
1 0.1
2 0.3


could be made using MATLAB's interp1 function:

Te = interp1([1,2],[0.1,0.3],1.1)


Sketch these measurements on a piece of paper, make the calculation of an estimate for T at x=1.1 by hand using only a calculator, and then compare your result with that from executing the above.

## 6.2. 1-D Interpolation

An estimate of speed at x = 1.9 given the following measurements

x speed
1 11
2 33


could be made using MATLAB's interp1 function:

Te = interp1([1,2],[11,33],1.1)


Sketch these measurements on a piece of paper, make the calculation of an estimate for speed at x=1.9 by hand using only a calculator, and then compare your result with that from executing the above.

## 6.3. 2-D Interpolation

Suppose that you have the following set of measurements of height and temperature at four grid points and want to estimate the height and temperature somewhere inside of the four points.

x   y   H   T
1   1   10  94
1   2   20  95
2   1   30  96
2   2   40  97


Write a program that uses interp2 to estimate the height at (x,y) = (1.5,1.5). Do the same for temperature.

## 6.4. 2-D Interpolation

Suppose that you have the following set of measurements of Bx and By at four grid points and want to estimate Bx and By somewhere inside of the four points.

x   y   Bz   Bz
1   1   1.1  0.4
1   2   2.2  0.5
2   1   3.3  0.6
2   2   4.4  0.7


Write a program that uses interp2 to estimate Bx and By at (x,y) = (1.35,1.35).

## 6.5. 2-D Interpolation - Special points

Suppose that you have a set of measurements of height at four grid points and want to estimate the height somewhere inside of the four points.

x   y   H
1   1   10
1   2   20
2   1   30
2   2   50


Note that the center point of (x,y) = (1.5,1.5) is special: the result of 2-D interpolation is the average of the diagonal averages. That is, the interpolated value at (x,y) = (1.5,1.5) is the average of all four corners (10+50+20+30)/4 = 27.5, which is (30+25)/2 = 27.5. At the center point, each of the corners contributes equal weight to the average.

Write a program that uses interp2 to estimate the value of H at (x,y) = (1.5,1.5) and verify that it gives 27.5. Next, modify the values in the H column and compare the result of using your program to the diagonal averaging method.

Points along the border of the square with corners given by the x and y values are also special in that 1-D interpolation can be used. Write a program that uses interp1 and interp2 along with a test point on the border to show that they both give the same result.

## 6.6. 2-D Interpolation

Suppose that you have the following set of measurements of height H and temperature T at nine grid points and want to estimate the height and temperature at a set of given points.

x   y   H  T
1   1   1  4
2   2   2  7
1   2   3  5
2   1   4  8
1   3   9  4
2   3   2  2
3   3   1  5
3   2   9  9
3   1   8  8



Draw a sketch of the grid points and their associated measurements.

Write a program that uses interp2 to estimate the height at

(x,y) = (1.75,1.25)

(x,y) = (3,1.5)

and

(x,y) = (3,1)

## 6.7. meshgrid

The function meshgrid can be used to create matrices that are required for the input to interp2.

For example, instead of writing

X  = [1    2;   1   2];
Y  = [1    1;   2   2];


one could write

[X,Y] = meshgrid([1:2],[1:2])


Use both methods to create a grid that corresponds to measurements made in the range of x=9-10 in increments of 0.5 and y=100-102 in increments of 0.5.

Both methods should give the same matrices X and Y.

## 6.8. Streamlines in 2-D

The values of velocity (in m/s) measured at the location where the measurements were made (in m) are shown in the following table.

 x   y   Vx  Vy
1   1   4   4
1   2   4   4
2   1   4   4
2   2   4   4


A massless particle is released from the point x=1.5 and y=1.5.

1. Draw the positions on a piece of paper and label the velocities.
2. Using a step size of 0.1 seconds, estimate the position of the particle after two steps.
3. Using a step size of 0.001 seconds, estimate the position of the particle after 100 steps.

## 6.9. Streamlines in 2-D

At the four points (x,y) = { (1,1), (2,1), (1,2), (2,2) }, wind velocity in the y direction is measured to be {1.0, 2.0, 3.0, 4.0}. The wind velocity in the x direction is {1.0, 1.0, 1.0, 1.0}. If a massless particle is placed at (x,y) = (1,1) and released, what will its position be a short time later (after a step of h=0.1)?

Sketch these measurements on a piece of paper and then make this calculation of the position after one step using only a calculator.

Calculate the position of the particle after two steps. Note that you may use MATLAB to do the interpolation calculation - just write out the code that you entered into MATLAB to do the calculation.

## 6.10. Streamlines in 2-D

At the four points (x,y) = { (1,1), (2,1), (1,2), (2,2) }, wind velocity in the y direction is measured to be {1.0, 1.0, 1.0, 1.0}. The wind velocity in the x direction is {1.1, 2.2, 3.3, 4.4}. If a massless particle is placed at (x,y) = (1,1) and released, what will its position be a short time later (after a step of h=0.01)?

Sketch these measurements on a piece of paper and then make this calculation of the position after one step using only a calculator.

Calculate the position of the particle after two steps. Note that you may use MATLAB to do the interpolation calculation - just write out the code that you entered into MATLAB to do the calculation.

## 6.11. Streamlines in 2-D

At the four points (x,y) = { (1,1), (2.5,1), (2.5,2.5), (1,2.5) }, wind velocity in the x direction is measured to be {1.0, 2.0, 3.0, 4.0}. The wind velocity in the y direction is {1.0, 1.0, 1.0, 1.0}. If a massless particle is placed at (x,y) = (1,1) and released, what will its position be a short time later (after a step of h=0.05)?

Notation: The first pair of values in the point list corresponds to the first value in the vx list and first value in the vy list. The second pair of values in the point list corresponds to the second value in the vx list and second value in the vy list, etc.

Sketch these measurements on a piece of paper and then make this calculation of the position after one step using only a calculator.

Calculate the position of the particle after two steps. Note that you may use MATLAB to do the interpolation calculation - just write out the code that you entered into MATLAB to do the calculation.

## 6.12. Streamlines in 2-D

At the four points (x,y) = { (1,1), (2.5,1), (2.5,2.5), (1,2.5) }, wind velocity in the x direction is measured to be {1.0, 2.0, 3.0, 4.0}. The wind velocity in the y direction is {1.0, 1.0, 1.0, 1.0}. If a massless particle is placed at (x,y) = (1,1) and released, what will its position be a short time later (after a step of h=0.05)?

Notation: The first pair of values in the point list corresponds to the first value in the vx list and first value in the vy list. The second pair of values in the point list corresponds to the second value in the vx list and second value in the vy list, etc.

Sketch these measurements on a piece of paper and then make this calculation of the position after one step using only a calculator.

Calculate the position of the particle after two steps. Note that you may use MATLAB to do the interpolation calculation - just write out the code that you entered into MATLAB to do the calculation.

## 6.13. Streamlines in 2-D

At the four points (x,y) = { (1,1), (2,1), (1,2), (2,2) }, wind velocity in the y direction is measured to be { 1.0, 2.0, 3.0, 4.0 }. The wind velocity in the x direction is {1.0, 1.0, 1.0, 1.0}. If a massless particle is placed at (x,y) = (1,1) and released, what will its position be a short time later (after 10 steps of h=0.01)? You may use the function interp2.

## 6.14. Streamlines in 2-D

Suppose that you have the following set of measurements of vx and vy at four grid points and want to estimate vx and vy somewhere inside of the four points.

x   y   vx   vy
1   1   1.1  0.4
1   2   2.2  0.5
2   1   3.3  0.6
2   2   4.4  0.7


Sketch out the measurements on a piece of paper and attempt to estimate vx and vy at (x,y) = (1.35,1.35).

Assume that vx and vy represent velocity measurements in m/s and a massless particle is released at (x,y) = (1.35,1.35). Estimate its position a short time, Δt = 0.01 seconds, later by making one step.

Save your program as HW5e.m.

Generalize your solution to the previous problem so that it computes the positions of the particle after making 100 steps, with each step of size Δt = 0.01. Note that you can plot estimated positions using

plot(px,py)


assuming px and py are arrays containing the positions of the particle for 50 steps (corresponding to 1 second of time).

## 6.15. Streamlines in 2-D

At the four points (x,y) = { (1,1), (2,1), (1,2), (2,2) }, wind velocity in the y direction is measured to be { 1.0, 2.0, 3.0, 4.0 }. The wind velocity in the x direction is {1.0, 1.0, 1.0, 1.0}. If a massless particle is placed at (x,y) = (1,1) and released, what will its position be a short time later (after 10 steps of h=0.01)? You may use the function interp2.

## 6.16. Streamlines in 2-D

First, compute and draw streamlines for the vector field created by executing

clear;
figure(1);clf
[X,Y] = meshgrid([-pi:0.1:pi],[-2:0.1:2]);
U = Y;
V = -sin(X);

quiver(X,Y,U,V);


Then determine the following:

1. For h=0.1, how many steps are required for the streamline to rotate 360o?
2. For h=0.01, how many steps are required for the streamline to rotate 360o?
3. Explain the difference in the path of the streamline for h=0.1 and h=0.01.

## 6.17. Streamlines in 2-D

The following program generates a hypothetical set of velocity measurements at locations on a 2-dimesional grid. The matrix U contains the velocities in the x-direction; the matrix V contains the velocities in the y-direction. The matrix X contains the x-positions of the measurements of velocity and Y contains the y-positions of the measurements of velocity.

clear;
figure(1);clf
[X,Y] = meshgrid([-pi:0.1:pi],[-2:0.1:2]);
U = Y;
V = -sin(X);


A massless particle is released from a location x,y = 1.0, 1.0.

### 6.17.1. Part I

Using a step size of h = 0.1, how many steps are required before the particle moves into a region where no measurements are available? What is the position of the particle just before it moves into this region? Answer these two questions also for h = 0.01.

Save your answer in a file named Midterm_4I.m. When I execute the program, I should see two plots, with each plot showing the trajectory of the particle. The title of each plot should list the exit position and the number of steps taken to get to the exit position.

### 6.17.2. Part II

Suppose that the velocities were measured on a finer grid. The following program creates a set of velocities with the same overall shape as the previous case, but there are now measurements every 0.05 units in x and 0.05 units in y (in Part I, the spacing was 0.1 units in x and 0.1 units in y).

clear;
figure(1);clf
[X,Y] = meshgrid([-pi:0.05:pi],[-2:0.05:2]);
U = Y;
V = -sin(X);


Using a step size of h = 0.1, how many steps are required before the particle moves into a region where no measurements are available? What is the position of the particle just before it moves into this region. Answer these two questions also for h = 0.01.

Save your answer in a file named Midterm_4II.m. When I execute the program, I should see two plots, with each plot showing the trajectory of the particle. The title of each plot should list the exit position and the number of steps taken to get to the exit position.

## 6.18. Streamlines in 2-D - Pendulum Interpretation

The following program generates a hypothetical set of velocity measurements at locations on a 2-dimensional grid. The matrix U contains the velocities in the x-direction; the matrix V contains the velocities in the y-direction. The matrix X contains the x-positions of the measurements of velocity and Y contains the y-positions of the measurements of velocity.

clear;
figure(1);clf
[X,Y] = meshgrid([-pi:0.1:pi],[-2:0.1:2]);
U = Y;
V = -sin(X);


Note that in this problem you are given the position and velocity matrices; in previous problems, you had to manually create these matrices.

You can visualize the vectors using quiver(X,Y,U,V).

A massless particle is released from a location (x,y) = (1.0, 1.0).

Using a step size of dt = 0.1, how many steps are required before the particle moves into a region where no measurements are available? What is the position of the particle just before it moves into this region? Answer these two questions also for dt = 0.01.

When I execute the program, I should see two plots, with each plot showing the trajectory of the particle (plot its y position versus its x position), one for dt = 0.1 and the other for dt = 0.01. The title of each plot should list the exit position and the number of steps taken to get to the exit position.

## 6.19. Streamlines in 2-D

Particles in a box with a random velocity field

You are given a file containing the positions where velocity measurements have been made. The velocity measurements have been made inside of a box with rigid walls. You are given a second file containing the initial positions of particles. When these particles are released, compute their trajectories.

### 6.19.1. Part I

Reading the files

Write a program that reads the data in the following two files into matrices. Do not use a text editor to modify the contents of the file prior to reading it into MATLAB.

1. velocities.txt - After your program is executed, you should have a matrix named Pos_and_Vel with four columns and 10201 rows.
2. positions.txt - After your program is executed, you should have a matrix named Particles with two columns and 9 rows.

### 6.19.2. Part II

Converting the input positions and velocities into a form that can be used by MATLAB's interp2 function

Similar to previous problems, you are given a table of x,y,vx,vy values. To use MATLAB's interp2 function, each of these variables must be represented as a matrix. However, after reading the files, you have only 1-D arrays containing these parameter values. For example,

x = Pos_and_Vel(:,1)  % Will have 101x101 = 10201 values.
y = Pos_and_Vel(:,2)  % Will have 101x101 = 10201 values.
vx = Pos_and_Vel(:,3) % Will have 101x101 = 10201 values.
vy = Pos_and_Vel(:,4) % Will have 101x101 = 10201 values.


are the position and velocity arrays. The next steps is to create matrices that correspond to this information. After running your code, you should have four matrices, X, Y, Vx, and Vy. Each of these matrices should have 101 rows and 101 columns. The position matrices should have the form

X

0 0.1 ... 10.0
0 0.1 ...
.
.
.
0 0.1 ... 10.0


Y

0.0 0.0 ... 0.0
0.1 0.1 ... 0.1
.
.
.
10 10 ... 10


And the Vx matrix should have the form

Vx

0 0 ... 0 0
0 a ... b 0
.
.
.
0 c ... d 0
0 0 ... 0 0


Where the value at

• a corresponds to the x-component of the velocity measured at x=0.1 and y=0.1,
• b corresponds to the x-component of the velocity measured at x=9.9 and y=0.1,
• c corresponds to the x-component of the velocity measured at x=0.1 and y=9.9, and
• d corresponds to the x-component of the velocity measured at x=9.9 and y=9.9.

### 6.19.3. Part III

Stepping the particles in time and visualizing their trajectories

Trace the path of a particle with an initial position of x=0.5 and y=0.5 until it steps outside of the grid or stops moving and then move on to part IV. Choose an appropriate step size Δt.

### 6.19.4. Part IV

Repeat part III for all of the particles with locations specified in the file positions.txt.

When executed, your program should show the traces of all particles on the same plot.

### 6.19.5. Part V

Modify your program so that vx = 1.0 and vy = 1.0 outside of the range $4.8 \le x \le 5.2$ and $4.8 \le y \le 5.2$.

Plot the trace of 10 particles that start at x = 5.0 and y = [4.5:0.1:5.4] until they step outside of the boundary or stop moving.

## 6.20. Dipole Field I

Charged particles in a static magnetic field appear to be guided by magnetic field lines. For example, in this video, the net motion of the charged particle is along the direction of the magnetic field. On average, the particle does not move in the horizontal direction.

As shown in the image, when iron filings are placed on a piece of paper near a magnetic field, the filings tend to align with the local magnetic field direction. From www.wired.com on April 14 2018 21:43:17.

These examples illustrate the usefulness of drawing streamlines associated with a magnetic vector field.

The program DipoleField.m creates magnetic field vectors associated with a dipole that is aligned in the z-direction. The magnetic field in the y-direction is zero.

Starting with the program DipoleField.m, write a program that draws a streamline that starts at the position (x,z) = (cos(π/4),sin(π/4)). In order to run the program, you also need to save B_dipole.m in the same directory as DipoleField.m

Extend this program so that on one plot streamlines are shown for starting positions of (x,z) = (cos(θ),sin(θ)) with θ = π/16, 2π/16, ..., 7π/16 (a total of 7 angles). Use the break command to stop the calculation of the streamline if it enters a position with x2 + y2 < 1 (or use a while loop).

The resulting image should look similar to that shown below. From raw.githubusercontent.com on May 19 2019 09:33:02.

## 6.21. Dipole Field II

Solve the previous problem using h = -0.1 (so that the streamlines go inwards). Modify h at each step so that its magnitude is inversely proportional to the magnitude of the length of the previous step.

## 6.22. Dipole Field III

Create a 3-D dipole field and plot outward going field lines outside of the sphere only for points that start on a unit sphere. Choose starting points of θ = π/16, 2π/16, ..., 7π/16 and φ = 0, π/16, 2π/16, ..., 15π/16.

## 6.23. Streamlines in 3-D

Create a 3-D vector field on your own and draw streamlines for it.