This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Engineering 101 Engineering 101
Lecture 25 12/4/07
3D Data and Mandelbrot Quote of the Day (Redux) Quote of the Day (Redux)
Those who will not economize will have to Confucius agonize. Project 7 – Revisiting time steps Project 7 – Revisiting time steps Suppose we have the following values: Now calculate some things: tFinal = 2, dt = 0.05, numPlots = 3, numSteps = 11 totalTimeSteps = ceil(tFinal/dt) = 40 plotSteps = ceil(totalTimeSteps/numPlots) = 14 2D Average Temp: 11 22 33 3D Temp Plot: 14, 28, 42 40 11 14 22 28 33 40 Note the increments for data logging: Note the final data logging profile (to use with multiSteps): Also, if the value for dt does not allow for a perfect match with tFinal, treat the step immediately after tFinal as the time value for the last plot tFinal = 2, dt = 0.03, numPlots = 3, numSteps = 11 t values for each step: 0, 0.03, 0.06, 0.09, …, 1.98, 2.01 (note 2.00 is missing) tempSim(2,50,0.05,8,11,3) versus tempSim(2,50,0.03,8,11,3) More Plotting in MATLAB More Plotting in MATLAB So far we have learned about normal xy plots You can actually plot multiple plots on the same axes using an alternative method hold on will prevent the graph from being redrawn and all subsequent plots will be plotted on top of old graphs hold off will turn off this feature Other TwoDimensional Plots Other TwoDimensional Plots
polar(theta, r) bar(x, y) barh(x, y) pie(x) A polar plot A vertical bar chart A horizontal bar chart (x = labels, y = values) A pie chart Three Dimensional Plots Three Dimensional Plots Let’s say you have some time dependent data For example the location of an object is described by the equations: x(t) = exp(0.2*t) * cos(2*t) y(t) = exp(0.2*t) * sin(2*t) We can plot this in 2D easily 1 Which MATLAB commands will plot the set of points visited by an object following the trajectory: x(t) = exp(0.2*t) * cos(2*t) 2 y(t) = exp(0.2*t) * sin(2*t) Exercise Exercise 3 4 1 Which MATLAB commands will plot the set of points visited by an object following the trajectory: x(t) = exp(0.2*t) * cos(2*t) 2 y(t) = exp(0.2*t) * sin(2*t) Exercise Exercise 3 4 plot2d.m Three Dimensional Plots Three Dimensional Plots
t = 0:0.1:10; x = exp(0.2*t) .* cos(2*t); y = exp(0.2*t) .* sin(2*t); plot(x,y); This produces a 2d plot, but it tells us nothing about time. plot3d.m Three Dimensional Plots Three Dimensional Plots
t = 0:0.1:10; x = exp(0.2*t) .* cos(2*t); y = exp(0.2*t) .* sin(2*t); plot3(x, y, t); Now we can immediately see the time dependence. Three Dimensional Plots Three Dimensional Plots 3D plots are particularly useful when we have data or functions that vary in more than one spatial dimension. For example let’s say we wanted to know what the function z = (x2 – y2) exp((x2 + y2)) looks like. Set Up a Mesh Set Up a Mesh Use the built in function: [x,y]= meshgrid(xstart:xinc:xend, ystart:yinc:yend) Once the grid is set up, then computing our function is simple: z = (x.^2 – y.^2) .* exp((x.^2 + y.^2)); Three Dimensional Plots Three Dimensional Plots The simplest way to plot is a pseudocolor plot pcolor(x, y, z); You can also plot contours contour(x, y, z); contourf(x, y, z); %filled contour Perspective views can be obtained with mesh(x, y, z); meshc(x, y, z); % with contour surf(x, y, z); surfc(x, y, z); % with contour surfl(x, y, z); % with lighting more3d.m Plotting in 3D Plotting in 3D
[x, y] = meshgrid(4:0.1:4, 4:0.1:4); z = (x.^2y.^2).*exp((x.^2+y.^2)); pcolor(x,y,z); or contour(x,y,z); mesh(x,y,z); surf(x, y, z); Options with 3D Plots Options with 3D Plots In shaded plots (e.g. pcolor, surf) you can choose different shading options: shading faceted shading flat shading interp % can be slow You can also choose different colormaps colormap pink colormap copper See the online help for more options Figure Window Tools Figure Window Tools The tools in the figure window can also be used to manipulate your image. These include a tool to zoom in/out, and to rotate and change the camera viewpoint. Making Movies in MATLAB Making Movies in MATLAB You can use the getframe and movie commands to create and play a movie. We will make a motion picture feature starring the thrilling function: cos ( t / 2 ) x − y cos( t ) + 2 xy sin ( t ) exp − x + y
2 2 2 2 {( ) } [( 2 )] [x, y] = meshgrid(4:0.1:4,4:0.1:4); f1 = exp((x.^2 + y.^2)); f2 = x.^2 y.^2; f3 = 2.*x.*y; Making Movies in MATLAB Making Movies in MATLAB [x, y] = meshgrid(4:0.1:4,4:0.1:4); f1 = exp((x.^2 + y.^2)); f2 = x.^2 y.^2; f3 = 2.*x.*y; Making Movies in MATLAB Making Movies in MATLAB
Rotate Twice for t = 0:pi/20:4*pi z = (cos(0.5*t)^2).*(f2.*cos(t)+f3.*sin(t)).*f1; end [x, y] = meshgrid(4:0.1:4,4:0.1:4); f1 = exp((x.^2 + y.^2)); f2 = x.^2 y.^2; f3 = 2.*x.*y; Making Movies in MATLAB Making Movies in MATLAB colormap copper; for t = 0:pi/20:4*pi z = (cos(0.5*t)^2).*(f2.*cos(t)+f3.*sin(t)).*f1; surfl(x,y,z); axis([4 4 4 4 0.5 0.5]); % axes always same shading interp; end surf_movie.m [x, y] = meshgrid(4:0.1:4,4:0.1:4); f1 = exp((x.^2 + y.^2)); f2 = x.^2 y.^2; f3 = 2.*x.*y; j=1; colormap copper; for t = 0:pi/20:4*pi z = (cos(0.5*t)^2).*(f2.*cos(t)+f3.*sin(t)).*f1; surfl(x,y,z); axis([4 4 4 4 0.5 0.5]); % axes always same shading interp; mov( j ) = getframe; j = j+1; end movie(mov); %play the movie Making Movies in MATLAB Making Movies in MATLAB Default Arguments in Functions Default Arguments in Functions You may have noticed that some of the bulitin MATLAB functions can accept and return different numbers of arguments. This is handled by using the values nargin and nargout. nargin/nargout nargin/nargout nargin is equal to the number of input arguments passed into the function. nargout is equal to the number of output arguments requested from the function. Using nargin Using nargin Let’s say I want to have a procedure that has some default arguments if none is specified. function val = simulate(time, startval) val = zeros(time); if nargin<2  isempty(startval) val(1) = 0; else val(1) = startval; end for t = 2:time val(t) = someFunction(val(t1)); end Using nargout Using nargout
function [m,v]=MeanVar(X) % MeanVar Computes the mean and variance n=size(X,1); m=mean(X); if nargout>1 temp=Xm(ones(n,1),:); v=sum(temp.*temp)/(n1); end 1 Which function will take in up to 3 matrices and return the average of these matrices? Exercise Exercise
2 3 4 1 Which function will take in up to 3 matrices and return the average of these matrices? Exercise Exercise
2 3 4 Programming Example: Programming Example: The Mandelbrot Set Programming Example: Programming Example: The Mandelbrot Set Fractals are of interest because, in addition to being mathematically beautiful objects, they have the property of selfsimilarity. In this respects they are like a number of natural and manmade systems like coastlines and rough surfaces. The Mandelbrot set The Mandelbrot set Generated from the equation: zm +1 = 2 zm +c The Mandelbrot set The Mandelbrot set Generated from the equation: zm +1 =
c=1 2 zm +c z1=0 z2=1 z3=2 z4=5 z5=2 z6=677 z7=458330 6 c=0.1 z1=0 z2=0.1 z3=0.11 z4=0.112 z5=0.11256 1 … The Mandelbrot set The Mandelbrot set We determine the value of a point by determining how many iterations it takes to grow greater than 2 and dividing by the total number of iterations. For example, consider 10 total iterations M(1) = 0.4 M(0.1) = 1 (it never becomes greater than 2) c=1 z1=0 z2=1 z3=2 z4=5 z5=2 z6=677 z7=458330 6 c=0.1 z1=0 z2=0.1 z3=0.11 z4=0.112 z5=0.11256 1 … The Mandelbrot set The Mandelbrot set The trick is we do this with complex numbers. So every point in the plane has a real value given by the x axis and an imaginary value given by the y axis. Squares of imaginary numbers are taken in the standard way (1+2i)2 = (1+2i)(1+2i) = 1+2i+2i+4i2 = 1+4i 4 = 3+4i The condition to stop the iteration will be that the norm of the number (the number times its complex conjugate) is greater than 2. The Mandelbrot set The Mandelbrot set So to make the set we must first write a function that will take as input: It will have to set the initial values of z for each c to 0. Then we must make a loop that will repeatedly apply the equation to each value of z a matrix c where each element of the matrix is a complex number a number of iterations, niters zm +1 = 2 zm +c The Mandelbrot set The Mandelbrot set
function res = mandelbrotIterate (c, niters) z = zeros(size(c)); for i = 1:niters z = z.^2 + c; end The Mandelbrot set The Mandelbrot set Since we only want to continue to iterate the z’s that are less than 2 we will use a logical array called active The value of active will be 1 if abs(z) is still less than or equal to 2 The Mandelbrot set The Mandelbrot set
function res = mandelbrotIterate (c, niters) z = zeros(size(c)); active = ones(size(c)); for i = 1:niters z(active) = z(active).^2 + c(active); active = abs(z) <= 2; end The Mandelbrot set The Mandelbrot set As a final step, we need to make sure than when a site first becomes inactive, the program determines the iteration in which it exceeded 2 and calculates the value of iteration/total iterations for the result, res. When a site first becomes inactive, the program must determine the When a site first becomes inactive, the program must determine the iteration in which it exceeded 2 and calculate the value of iteration/total iterations for the result, res. Which function does this? 1 2 When a site first becomes inactive, the program must determine the When a site first becomes inactive, the program must determine the iteration in which it exceeded 2 and calculate the value of iteration/total iterations for the result, res. Which function does this? 1 2 The Mandelbrot Set The Mandelbrot Set
function res = mandelbrotIterate (c, niters) z = zeros(size(c)); active = ones(size(c)); res = ones(size(c)); for i = 1:niters z(active) = z(active).^2 + c(active); newactive = abs(z) <= 2; res(active & ~newactive)=i/niters; active = newactive; end The Mandelbrot Set The Mandelbrot Set Now we have a function mandelbrotIterate(c, niters) that will calculate the result for a matrix of complex values c over a number of iterations niter Next we must design a function that will create the matrix c, call mandelbrotIterate and plot the result in a pseudocolor plot. We will call this function drawMandelbrot and pass it the range [xmin xmax ymin ymax], the resolution (number of boxes in each direction) and the number of iterations (niter). The Mandelbrot Set The Mandelbrot Set
function drawMandelbrot(range,resolution,niter) % we will have to create c here res = mandelbrotIterate(c,niter); pcolor(res); shading flat; colormap jet; Which function fragment creates c properly? Recall: Which function fragment creates c properly? Recall: 1 2 3 4 range is a vector that holds [xmin xmax ymin ymax] resolution is the number of boxes in each direction Each box in c holds the complex number x+iy Which function fragment creates c properly? Recall: Which function fragment creates c properly? Recall: 1 2 3 4 range is a vector that holds [xmin xmax ymin ymax] resolution is the number of boxes in each direction Each box in c holds the complex number x+iy The Mandelbrot Set The Mandelbrot Set
function drawMandelbrot(range,resolution,niter) xvals = range(1):(range(2)range(1))/(resolution1):range(2); yvals = range(3):(range(4)range(3))/(resolution1):range(4); [re, im] = meshgrid(xvals, yvals); c = re + 1.0i.*im; res = mandelbrotIterate(c,niter); pcolor(xvals, yvals, res); shading flat; colormap jet; Now if we execute the command drawMandelbrot([2.5 2.5 2 2],400,40) we obtain the image: Mandelbrot set Mandelbrot set Code Available on CTools Code Available on CTools Code for mandelbrotIterate, drawMandelbrot and mandelbrotgui are available on the course website resources in the MATLAB folder. Now for some cool pictures… Now for some cool pictures… Since MATLAB’s resolution is so bad, let’s have the wonderful contributors at wikipedia makes the plots for us! And, of course, we need some music 1 : 60,000,000,000 Ordinary Monitor : 20,000,000 KM Next Lecture Next Lecture MATLAB I/O and a bit more on GUIS ...
View
Full
Document
This note was uploaded on 03/12/2010 for the course ENGIN 101 taught by Professor Jeffringenberg during the Fall '07 term at University of Michigan.
 Fall '07
 JeffRingenberg

Click to edit the document details