Streamlines
From ComputingForScientists
1. References
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 [2]. Scalars are typically visualized using contour lines and images with colorbars, both of which are shown in [3].
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 [4] 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 ydirections is proportional to the velocities in the x and ydirection, 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 zoomedin 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, bilinear 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. 1D
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 inbetween. 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 1D 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 1D interpolation was needed to estimate vy, because vy only varied in the xdirection.
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. 2D
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 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 rearrange the order of the values in X
and Y
provided that we reordered the values in H
for consistency:
X = [ 1 , 1 ; 2 , 2] Y = [ 1 , 2 ; 1 , 2] H = [10 , 20 ; 30 , 40]
The above setup 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. 2D and 1D equivalence
Sometimes a problem that appears to require 2D interpolation can actually be solved using 1D 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 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 ypositions 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(t1). Last velocity is vx(t1). px(t) = px(t1) + dt*vx(t1); % Take a step in y. New position is py(t). Last position is py(t1). Last velocity is vy(t1). py(t) = py(t1) + dt*vy(t1); % 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
Thus far, we have considered cases where F_{x} and F_{y} 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 = F_{x}(x,y)
(2) dy / dt = F_{y}(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 that is made by stepping a distance dx in the xdirection and dy in the ydirection. This triangle is similar to a triangle formed with F_{x} in the xdirection, F_{y} in the ydirection, and F as the hypotenuse so that
dx / ds = F_{x}(x,y) /  F 
dy / ds = F_{y}(x,y) /  F 
To determine the x and y values that correspond to a line that is always in the direction of , 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 and 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 and so that
Δx / Δs = F_{x}(x,y) /  F 
Δy / Δs = F_{y}(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 , 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. 1D 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.
Answer 

clear % Hand calculation using similar triangle method. % Te = 0.1 + (0.2)*(x1.0)/(2.01.0) % plug in x = 1.1 Te = 0.1 + (0.2)*(1.11.0)/(2.01.0) % Hand calculation by solving for slope and intercept: % First, find m and b in % T = mx+b. % (1) 0.1 = m*1.0 + b % (2) 0.3 = m*2.0 + b % Solve (1) for b % (3) b = 0.1m*1.0 % Plug (3) into (2) % (4) 0.3 = m*2.0 + 0.1m*1.0 % Solve for m % (5) m = 0.2/(2.01.0) % plug (5) into (3) % (6) b = 0.1  (0.2)*1.0/(2.01.0) % (7) b = 0.1 (verify by drawing line and checking intercept) % % Equation of line is % T = 0.2*x  0.1 % when x = 1.1, % Te = 0.2*1.1  0.1 = 0.12 figure(1);clf plot([1,2],[0.1,0.3],'k.','MarkerSize',20); hold on; axis([0 3 0 0.5]) xlabel('x') ylabel(gca, 'y','Rotation',0); set(gca,'XTick',[1,2]) grid on; Te = interp1([1,2],[0.1,0.3],1.1) plot(1.1,Te,'r.','MarkerSize',20) text(1.2,Te,sprintf('Estimate of T at x=1.1 is T_e = %.2f',Te)) 
6.2. 1D 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. 2D 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.
Answer 

% Important: Sketch out values so that you can check that % value from interp2 makes sense! X = [1 2;1 2]; Y = [1 1;2 2]; H = [10,30;20,40]; T = [94,96;95,97]; % Diagonal averages are 15 and 25. Average of diagonals is 20. % The following should give 20. He = interp2(X,Y,H,1.5,1.5) % Diagonal averages are 95.5 and 95.5. Average of diagonals is 95.5. % The following should give 95.5. Te = interp2(X,Y,T,1.5,1.5) figure(1);clf plot([1,1,2,2],[1,2,1,2],'k.','MarkerSize',20) grid on; hold on; plot(1.5,1.5,'x') text(1.55,1.5,sprintf('H_e = %.2f',He)) text(1.2,1,'H = 10') text(1.2,2,'H = 20') text(2.2,1,'H = 10') text(2.2,2,'H = 40') axis([0 3 0 3]) xlabel('x') ylabel(gca, 'y','Rotation',0); set(gca,'XTick',[1,2]) set(gca,'YTick',[1,2]) axis square figure(2);clf plot([1,1,2,2],[1,2,1,2],'k.','MarkerSize',20) grid on; hold on; plot(1.5,1.5,'x') text(1.55,1.5,sprintf('T_e = %.2f',Te)) text(1.2,1,'H = 94') text(1.2,2,'H = 95') text(2.2,1,'H = 96') text(2.2,2,'H = 97') axis([0 3 0 3]) xlabel('x') ylabel(gca, 'y','Rotation',0); set(gca,'XTick',[1,2]) set(gca,'YTick',[1,2]) axis square X = [1 2;1 2]; Y = [1 1;2 2]; H = [10,10;20,40]; interp2(X,Y,H,1.5,1.5) 
6.4. 2D 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. 2D 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 2D 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 1D 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.
Answer 

X = [1 2;1 2]; Y = [1 1;2 2]; H = [10,30;20,50] He = interp2(X,Y,H,1.5,1.5) % Test value along y = 1.0. He2 = interp2(X,Y,H,1.5,1.0) % Here the problem reduces to 1D interpolation. % The x values are the same, but the y values % passed to interp1 are the H values. He2 = interp1([1,2],[10,30],1.5) figure(1);clf plot([1,1,2,2],[1,2,1,2],'k.','MarkerSize',20) grid on; hold on; plot(1.5,1.5,'r.','MarkerSize',20) plot(1.5,1.0,'g.','MarkerSize',20) text(1.45,1.6,sprintf('H_e = %.2f',He),'Color','r') text(1.45,1.1,sprintf('H_e = %.2f',He2),'Color','g') text(1.1,1,'H = 10') text(1.1,2,'H = 20') text(2.1,1,'H = 30') text(2.1,2,'H = 50') axis([0 3 0 3]) xlabel('x') ylabel(gca, 'y','Rotation',0); set(gca,'XTick',[1,2]) set(gca,'YTick',[1,2]) axis square % Test value along x = 1.0. He3 = interp2(X,Y,H,1.0,1.3) He3 = interp1([1,2],[10,20],1.3) 
6.6. 2D 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)
Answer 

clear [X,Y] = meshgrid([1:3],[1:3]); H = [1 4 8;3 2 9;9 2 1]; T = [4 8 8;5 7 9;4 2 5]; H1 = interp2(X,Y,H,1.75,1.25) T1 = interp2(X,Y,T,1.75,1.25) H2 = interp2(X,Y,H,3,1.5) T2 = interp2(X,Y,T,3,1.5) H3 = interp2(X,Y,H,3,1) T3 = interp2(X,Y,T,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=910 in increments of 0.5 and y=100102 in increments of 0.5.
Both methods should give the same matrices X
and Y
.
Answer 

clear X1 = [9 9.5 10; 9 9.5 10; 9 9.5 10; 9 9.5 10; 9 9.5 10] Y1 = [100 100 100; 100.5 100.5 100.5; 101 101 101; 101.5 101.5 101.5; 102 102 102] [X,Y] = meshgrid([9:0.5:10],[100:0.5:102]) 
6.8. Streamlines in 2D
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
.
 Draw the positions on a piece of paper and label the velocities.
 Using a step size of
0.1
seconds, estimate the position of the particle after two steps.  Using a step size of
0.001
seconds, estimate the position of the particle after 100 steps.
Answer 

%X = [1 2;1 2]; % Don't need because no interpolation needed. %Y = [1 1;2 2]; x(1) = 1.5; y(1) = 1.5; h = 0.1; x(2) = x(1) + h*4; y(2) = y(1) + h*4; x(3) = x(2) + h*4; y(3) = y(2) + h*4; % etc. Substitute x(2) into equation for x(3) to get x(3) = x(1) + h*4 + h*4; Substitute x(3) into equation for x(4) to get x(4) = x(1) + h*4 + h*4 + h*4; % etc. In general, for Nsteps x(Nsteps+1) = x(1) + h*4*Nsteps; y(Nsteps+1) = y(1) + h*4*Nsteps; % So after 100 steps x(101) = x(1) + h*4*100; % 1.5 + 0.4 y(101) = y(1) + h*4*100; % 1.5 + 0.4 
6.9. Streamlines in 2D
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.
Answer 

px(2) = 1.10 py(2) = 1.10 vx(2) = 1.00 vy(2) = 1.30 px(3) = 1.20 py(3) = 1.23 vx(3) = 1.00 vy(3) = 1.66 X = [1 , 2 ; 1 , 2] Y = [1 , 1 ; 2 , 2] U = [1.0, 1.0; 1.0, 1.0] V = [ 1.0, 2.0; 3.0, 4.0 ] px(1) = 1.0; py(1) = 1.0; vx(1) = interp2(X,Y,U,px(1),py(1)); vy(1) = interp2(X,Y,V,px(1),py(1)); h = 0.1; Nsteps = 3; for t = 2:Nsteps px(t) = px(t1) + h*vx(t1); py(t) = py(t1) + h*vy(t1); % 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 
6.10. Streamlines in 2D
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.
Answer 

h = 0.01; X = [1 2 ;1 2 ]; Y = [1 1 ;2 2 ]; Vx = [1.1 2.2;3.3 4.4]; Vy = [1.0 1.0;1.0 1.0]; px(1) = 1.0; py(1) = 1.0; vx(1) = 1.1; vy(1) = 1.0; px(2) = px(1) + h*vx(1); py(2) = py(1) + h*vy(1); vx(2) = interp2(X,Y,Vx,px(2),py(2)); vy(2) = interp2(X,Y,Vy,px(2),py(2)); px(3) = px(2) + h*vx(2); py(3) = py(2) + h*vy(2); px py 
6.11. Streamlines in 2D
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.
Answer 

px(2) = 1.05 py(2) = 1.05 vx(2) = 1.10 vy(2) = 1.00 px(3) = 1.10 py(3) = 1.10 vx(3) = 1.20 vy(3) = 1.00 X = [1 , 2.5 ; 1 , 2.5] Y = [1 , 1 ; 2.5 , 2.5] U = [ 1.0, 2.0; 3.0, 4.0 ] V = [1.0, 1.0; 1.0, 1.0] px(1) = 1.0; py(1) = 1.0; vx(1) = interp2(X,Y,U,px(1),py(1)); vy(1) = interp2(X,Y,V,px(1),py(1)); h = 0.05; Nsteps = 3; for t = 2:Nsteps px(t) = px(t1) + h*vx(t1); py(t) = py(t1) + h*vy(t1); % 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 
6.12. Streamlines in 2D
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 2D
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
.
Answer 

% 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 ] px(1) = 1.0; py(1) = 1.0; vx(1) = interp2(X,Y,U,px(1),py(1)); vy(1) = interp2(X,Y,V,px(1),py(1)); % If px(1) and py(1) are 1.0, interpolation is not needed % because these positions are at a measurement point. % We could just say 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 % other points without having to manually change vx(1) and vy(1). h = 0.01; Nsteps = 10; for t = 2:Nsteps % Take a step. New position is px(t). Last position is px(t1). Last velocity is vx(t1). px(t) = px(t1) + h*vx(t1); py(t) = py(t1) + h*vy(t1); % 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)); 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 
6.14. Streamlines in 2D
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
.
Answer 

x = [1,2;1,2] y = [1,1;2,2] vx = [1.1,3.3;2.2,4.4] vy = [0.4,0.6;0.5,0.7] px(1) = 1.35; py(1) = 1.35; % Compute estimate of velocity at initial position vxi(1) = interp2(x,y,vx,px(1),py(1)) vyi(1) = interp2(x,y,vy,px(1),py(1)) Dt = 0.01; % New x,y positions is are position + estimated x,y velocities * Dt. px(2) = px(1) + Dt*vxi(1); py(2) = py(1) + Dt*vyi(1); 
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).
Answer 

Note that that px and py are NaN after step 24 because the particle has stepped outside of the region where measurements are available. clear x = [1,2;1,2] y = [1,1;2,2] vx = [1.1,3.3;2.2,4.4] vy = [0.4,0.6;0.5,0.7] px(1) = 1.35; py(1) = 1.35; Dt = 0.01; % New x,y positions is are position + estimated x,y velocities * Dt. for i = 1:99 % Compute estimate of velocity at initial position vxi(i) = interp2(x,y,vx,px(i),py(i)); vyi(i) = interp2(x,y,vy,px(i),py(i)); px(i+1) = px(i) + Dt*vxi(i); py(i+1) = py(i) + Dt*vyi(i); end plot(px,py,'.','MarkerSize',20) xlabel('x'); ylabel('y'); grid on; 
6.15. Streamlines in 2D
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 2D
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:
 For
h=0.1
, how many steps are required for the streamline to rotate 360^{o}?  For
h=0.01
, how many steps are required for the streamline to rotate 360^{o}?  Explain the difference in the path of the streamline for
h=0.1
andh=0.01
.
Answer 

1. and 2. 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); px(1) = 1.0; py(1) = 1.0; % Find interpolated velocities at initial point. vx(1) = interp2(X,Y,U,px(1),py(1)); vy(1) = interp2(X,Y,V,px(1),py(1)); h = 0.1; Nsteps = 77; % Value was guessed by trialand error for t = 2:Nsteps % Take a step px(t) = px(t1) + h*vx(t1); py(t) = py(t1) + h*vy(t1); % Find interpolated velocities at new position vx(t) = interp2(X,Y,U,px(t),py(t)); vy(t) = interp2(X,Y,V,px(t),py(t)); 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 ;hold on; plot(px,py,'g','LineWidth',3) plot(px,py,'k.','MarkerSize',10); plot(px(1),py(1),'kx','MarkerSize',15); 4. The actual path should close in on itself. With a smaller step size, the computed trajectory is more accurate (in reality, the algorithm used for stepping is poor. Most programs that compute streamlines use a more sophisticated stepping algorithm (e.g., RungeKutta) that than this one, which is called Forward Euler). 
6.17. Streamlines in 2D
The following program generates a hypothetical set of velocity measurements at locations on a 2dimesional grid. The matrix U
contains the velocities in the xdirection; the matrix V
contains the velocities in the ydirection. The matrix X
contains the xpositions of the measurements of velocity and Y
contains the ypositions 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.
Answer 


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.
Answer 


6.18. Streamlines in 2D  Pendulum Interpretation
The following program generates a hypothetical set of velocity measurements at locations on a 2dimensional grid. The matrix U
contains the velocities in the xdirection; the matrix V
contains the velocities in the ydirection. The matrix X
contains the xpositions of the measurements of velocity and Y
contains the ypositions 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.
Answer 

Many students wondered why reducing There is a mathematical way of showing that the trajectories close in on themselves. The vector field for this problem U = Y; V = sin(X); is related to the pendulum problem, which has equations of motion of dθ / dt = ω dω / dt = − sin(θ) If we divide the two equations to eliminate
If we let dθ be a vector in the xdirection and Failed to parse (unknown error\d): \d\omega be a vector in the y − direction, then we can write
and
If these vectors are plotted at various positions, they represent the direction the particle moves in the socalled "phase space" at those various positions. Consider the quiver plot. If the particle is released at 0, 0, it will never move. This is equivalent to starting the pendulum with an initial position of zero and an initial velocity of zero. Because the vector field given in this problem is related to the pendulum problem, we know that the trajectories should close on themselves in phase space. If the pendulum is started with a given initial θ and ω, it should eventually have those same values at some point in time because the motion of the pendulum is periodic and there is no damping. % First, compute and draw streamlines for the vector field created by executing [X,Y] = meshgrid([pi:0.1:pi],[2:0.1:2]); U = Y; V = sin(X); [X,Y] = meshgrid([pi:0.1:pi],[2:0.1:2]); U = Y; V = sin(X); px(1) = 1.0; py(1) = 1.0; % Find interpolated velocities at initial point. vx(1) = interp2(X,Y,U,px(1),py(1)); vy(1) = interp2(X,Y,V,px(1),py(1)); h = 0.1; t = 2; while (~isnan(px(t1)) && ~isnan(py(t1))) % Take a step px(t) = px(t1) + h*vx(t1); py(t) = py(t1) + h*vy(t1); % Find interpolated velocities at new position vx(t) = interp2(X,Y,U,px(t),py(t)); vy(t) = interp2(X,Y,V,px(t),py(t)); 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)); t = t + 1; end figure(1);clf; quiver(X,Y,U,V) hold on; plot(px,py,'g','LineWidth',3) plot(px,py,'k.','MarkerSize',10); plot(px(1),py(1),'kx','MarkerSize',15); % t2 for Nsteps because if px(2) is NaN (exit position) then one step % was taken. But t is incremented after exit position, so need % to subtract another one from t. ts = sprintf('\\Delta t = %.2f; N_{steps} = %d, (x,y)=(%.1f,%.1f)',... h,t1,px(t2),py(t2)); title(ts) h = 0.01; % Note  everything below is a repeat from above. % It would be better programming practice to do a % for loop for two values of h. To make the solution % easier to read, this was not done. t = 2; while (~isnan(px(t1)) && ~isnan(py(t1))) % Take a step px(t) = px(t1) + h*vx(t1); py(t) = py(t1) + h*vy(t1); % Find interpolated velocities at new position vx(t) = interp2(X,Y,U,px(t),py(t)); vy(t) = interp2(X,Y,V,px(t),py(t)); 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)); t = t + 1; end figure(2);clf; quiver(X,Y,U,V); hold on; plot(px,py,'g','LineWidth',3) plot(px,py,'k.','MarkerSize',10); plot(px(1),py(1),'kx','MarkerSize',15); ts = sprintf('\\Delta t = %.2f; N_{steps} = %d, (x,y)=(%.1f,%.1f)',... h,t1,px(t2),py(t2)); title(ts) 
6.19. Streamlines in 2D
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.
 velocities.txt  After your program is executed, you should have a matrix named
Pos_and_Vel
with four columns and 10201 rows.  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 1D 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 xcomponent of the velocity measured atx=0.1
andy=0.1
, 
b
corresponds to the xcomponent of the velocity measured atx=9.9
andy=0.1
, 
c
corresponds to the xcomponent of the velocity measured atx=0.1
andy=9.9
, and 
d
corresponds to the xcomponent of the velocity measured atx=9.9
andy=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 and .
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.
Answer 

clear fid = fopen('velocities.txt','r'); for i = 1:3 line = fgetl(fid); end for i = 1:101*101 line = fgetl(fid); if ~ischar(line); break end tmp = strrep(line,'(bottom boundary)','0'); tmp = strrep(tmp,'(left boundary)','1'); tmp = strrep(tmp,'(inner point)','2'); tmp = strrep(tmp,'(right boundary)','3'); tmp = strrep(tmp,'(top boundary)','4'); tmp = tmp(1:end); % or % tmp = line(1:36); % Pos_and_Vel(i,:) = str2num(tmp); Pos_and_Vel(i,:) = str2num(tmp); end fclose(fid); fid = fopen('positions.txt','r'); for i = 1:2 line = fgetl(fid); end for i = 1:9 line = fgetl(fid); if ~ischar(line); break end Positions(i,:) = str2num(line); end fclose(fid); k = 1; for i = 1:101 for j = 1:101 Vx(i,j) = Pos_and_Vel(k,3); Vy(i,j) = Pos_and_Vel(k,4); k = k + 1; end end [X,Y] = meshgrid([0:0.1:10],[0:0.1:10]); dt = 0.01; clf;hold on; for j =1:9 x(1) = Positions(j,1); y(1) = Positions(j,2); for i = 1:100 vxi = interp2(X,Y,Vx,x(i),y(i)); vyi = interp2(X,Y,Vy,x(i),y(i)); x(i+1) = x(i)+dt*(vxi); y(i+1) = y(i) +dt*(vyi); plot(x,y,'.','LineWidth',3); end end % Part V for i = 1:101 for j = 1:101 if x(k) > 4.8 && x(k) < 5.2 && y(k) > 4.8 && y(k) < 5.2 Vx(i,j) = vx(k)+10; Vy(i,j) = vx(k)+10; else Vx(i,j) = 1; Vy(i,j) = 1; end k = k + 1; end end 
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.
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 zdirection. The magnetic field in the ydirection 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 x^{2} + y^{2} < 1 (or use a while
loop).
The resulting image should look similar to that shown below.
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 3D 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 3D
Create a 3D vector field on your own and draw streamlines for it.