Unformatted text preview: University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz E7, Assignment 12
SOLUTION KEY Assigned: Tuesday, December 3, 2008 Due: 12:00 pm (noon), Wednesday, December 10, 2008 This assignment is on statistics, numerical differentiation, integration, and solution of ordinary differential equations. As before, turn in the hard copy of your published file to the drop boxes in Etcheverry 1109 and upload the soft copies of your script and your functions (the Mfiles) to bspace. Do not forget to name your main Mfile as lastname_firstname_SID_lab12.m NOTE: Do not forget to display the contents of your userdefined functions using the command type. MATLAB commands* introduced in this assignment: mean, std, bar, erf, erfc, quad 1. A package delivery company (such as FeDex or UPS) has 500 planes departing from one of its hubs and each plane is loaded with 10,000 packages. The weight of each package can be considered as a uniformly distributed independent random variable between 2 kg and 10 kg. a) Calculate the mean and standard deviation of the weight of each package. Hint: Define the weight of a package by W. Since each weight W is uniformly distributed between 2 kg and 10 kg, then the probability distribution function (PDF) of the package weigh is σ 2W = ∞ −∞ ∫ (w − μ W ) 2 pW ( w)dw = ∫ w2 pW ( w)dw − μ 2W =
−∞ ∞ 1 w2 dw − 62 = ? (10 − 2) ∫ 2 10 Complete the right hand side of the second equation.
>> muW_I=(10+2)/2 >> vW_I= 1/8*(10^3/32^3/3)6^2 >> sigmaW_I=sqrt(vW_I) b) Calculate the mean and standard deviation of the total weight L that each plane must carry. Hint: Since each plane carries 10,000 packages, and the weight of each package is an independent random variable with mean μW and standard deviation σ W , then, as discussed in Lecture Notes 23  Statistics, the mean and standard deviation of the total plane load L, μ L and σ L , are μ L = 10000 ⋅ μW and σ L = 10000 ⋅ σ 2W >> muL_I=10000*muW_I >> sigmaL_I=sqrt( 10000*sigmaW_I^2) c) The following MATLAB command generates a 500×1 array L, where the kth element
L(k) is a randomly generated number that simulates the load Lk carried by the kth plane.
⎧1 ⎪ for 2 ≤ w ≤ 10 pw ( w) = ⎨ 8 ⎪ 0 elsewhere ⎩
>> L = sum(8*rand(500,10000)+2,2); Therefore, as discussed in Lecture Notes 23  Statistics, the mean μW and standard deviation σ W of each package are respectively given by i) Use the MATLAB commands mean and std to respectively determine the sample mean and standard deviations
L= 1 500 ∑ L(k ) 500 k =1 and σ L = 1 500 ∑ ( L( k ) − L ) 2 500 k =1 μW = ∞ −∞ ∫ wp W ( w)dw = 1 1 102 − 22 10 + 2 ∫ wdw = 2 10 − 2 = 2 = 6 (10 − 2) 2 10 and respectively compare them with theoretical values μ L and σ L calculated above.
>> L = sum(8*rand(500,10000)+2,2); * Please refer to MATLAB help to learn how to use the functions introduced in this assignment. >> muL = mean(L) E7 Assignment 12 1 Assignment 12 2 E7 University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz >> sigmaL = std(L) 2. a) Write a function called myTrapIntegrator that computes the integral of a function using the trapezoidal method. This function should have 4 input parameters: • • • Function handle of the integrand function (funcHand) Lower limit of integration (a) Upper limit of integration (b) Number of subinterval endpoints used in the interval ( a ≤ x ≤ b ) for the trapezoidal method (N). The function should return the integral value for the input function with limits a, and b computed with the trapezoidal method with N1 subintervals. With N = 101 points, test your function by evaluating the following integral
π /2 ii) Use the function my_PDF shown below
function [PDF, x] = my_PDF(X,x) [z,x]=hist(X,x); PDF = (z/sum(z))*(length(x)/(max(X)min(X))); to determine the approximate PDF of L, as described in Lecture Notes 23  Statistics (start with 20 bins), and plot it using the bar command. iii) On the same figure plot a Gaussian PDF with mean μ L and standard deviation σ L , and comment on how well does the approximate PDF of L approximate a Gaussian distribution.
>> [pdf, x] = my_pdf(L,20); >> figure >> bar(x,pdf) >> hold on >> muX = muL_I; >> stdX = sigmaL_I; >> gaussian=@(x) (1/stdX/sqrt(2*pi))*... exp(((xmuX)/stdX).^2/2); • I1 = ∫ cos( x) + x
0 1/ 2 dx Use format long to display the result obtained with a high degree of precision. Solution: function [out] = myTrapIntegrator(funcHand,a,b,N) x = linspace(a,b,N); >> plot(x,gaussian(x),'go') >> hold off d) Assume now that the load carried by each plane is a Gaussian random variable with mean μ L and standard deviation σ L . As described in Lecture Notes 23 – Statistics and chapter 6.2 of the text, use the MATLAB erf and/or erfc functions to estimate the probability that a plane will carry a load larger than 60,400 Kg.
>> pR =1/2*erfc((60400muL_I)/sqrt(2)/sigmaL_I) dx = (ba)/(N1); out = 0; for i = 1:N1 % Number of intervals is (N1) out = out + .5*dx*(funcHand(x(i)) + funcHand(x(i+1))); end >> format long Assignment 12 3 E7 Assignment 12 4 E7 University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz >> f1 = @(x) cos(x) + x.^(1/2); >> myTrapIntegrator(f1,0,pi/2,101) >> It = myTrapIntegrator(f1,2,2,501) d) Until now, you have written integration functions that take the number of subinterval endpoints as an input to accurately compute the integral. In certain cases, you may not know a priori the number of subintervals required to accurately compute the integral. An approach to dealing with such cases is to increment the number of subintervals, and continue calculating the integral until the difference between the successive values of the integral is less than a prescribed tolerance. Write a function myTrapIntTol that uses this procedure to compute a given integral to within a given tolerance using the function myTrapIntegrator you wrote in Part (a). The function myTrapIntTol should take the following inputs: • • • • Function handle of the integrand function. (funcHand) Lower limit of integration (a) Upper limit of integration (b) Tolerance of Integration (tol) b) Calculate the relative error, ε(N), for the following set of N values, N = 6, 11, 21, 51, 101, 151, 251, 501, 1001 points The relative error is given by: ε (N ) = I − I1 ( N ) I where I is the exact value of the integral, which in this case can be calculated by hand, and I1 ( N ) is the value of the integral calculated using the trapezoidal method with a given values of N. Store the values of the relative errors obtained for all the values of N tested in a double array named trap.
>> exact = 1 + (pi/2)^1.5/1.5; >> trap = ; >> n = [6 11 21 51 101 151 251 501 1001] >> for N = n; trap=[trap(abs(myTrapIntegrator(f1,0,pi/2,N)exact)... /exact)]; end The function should start with 10 subintervals (11 subinterval endpoints) and increase the number of intervals by a factor of 2 (i.e., to 20, 40, 80, etc), until the computed value of the integral found is within a tolerance of tol from the value of the integral computed with half of the subintervals. Hint: use the definition of ε(N)
I1 ( N ) − I1 ( N 2 ) I1 ( N ) c) Evaluate the integrals I 2 using the quad function in MATLAB. I 2 = ∫ e − x dx
2 2 ε (N) = −2 Make sure that the function handle that is sent to the quad function is vectorized. Use
format long to display a higher degree of precision in your answers. In addition, Calculate and display the answers to the following integrals to the specified tolerances. Remember to use format long to show your answer to a high degree of precision. I 2 = ∫ e − x dx
2 compute the integral using myTrapIntegrator you created in part (a) with 501 subinterval endpoints. Display the answer obtained from each method.
>> format long >> f1 = @(x) exp(x.^2); >> Iq = quad(f1,2,2) 2 to tolerances of 105, and 1012 to tolerances of 102, and 1010 −2 π I3 = −π ∫ sin 2 ( x) + ecos( x ) dx Assignment 12 5 E7 Assignment 12 6 E7 University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz Solution:
end function [out] = myTrapIntTol(funcHand, a, b, tol) yy(i1)=(y1(i)y1(i1))/(x(i)x(i1)); >> yy2=cos(xx); >> subplot(3,1,2), plot(xx,yy,xx,yy2) N = 10; I = myTrapIntegrator(funcHand,a,b,N+1); while abs((ImyTrapIntegrator(funcHand,a,b,N/2+1))/I) > tol N = 2*N; I = myTrapIntegrator(funcHand,a,b,N+1); end out = I; >> title('First derivative of Sine') >> legend('sin','cos') >> for j=2:length(xx) p(j1)=xx(j1)+0.5*(xx(j)xx(j1)); q(j1)=(yy(j)yy(j1))/(xx(j)xx(j1)); end >> qq = cos(p); >> subplot(3,1,3), plot(p,q,p,qq) >> format long >> f1 = @(x) exp(x.^2); >> f2 = @(x) exp(cos(x)) +(sin(x)).^2; >> myTrapIntTol(f1,2,2,1e5) >> myTrapIntTol(f2,pi,pi,1e2) >> title('Second derivative of Sine') >> legend('sin','cos') 4. Load the rocket.mat file available on bspace to your workspace. The variables Time and Altitude will be loaded. These arrays are data collected from a twostage rocket. a) Plot the altitude vs. time
>> load 'rocket.mat' >> plot(Time,Altitude) >> title('Altitude vs. Time'),xlabel('Time'), ylabel('Altitude') 3. On a 3×1 array of subplots display the Matlab function sin and its first two derivatives, computed using finite differences. Compute the derivatives at 20 points between 0 and 2π. On each of the derivatives subplot, also plot the Matlab function cos. Clearly label your figure using legends.
>> x=linspace(0,2*pi,22); >> y1=sin(x); >> y2=cos(x); >> subplot(3,1,1), plot(x,y1,x,y2) >> title('Sine') >> legend('sin','cos') >> for i=2:length(x) xx(i1)=x(i1)+0.5*(x(i)x(i1)); b) Use finitedifference differentiation to compute the velocity and plot it.
>> for i=2:length(Time) Velocity(i1)=(Altitude(i)Altitude(i1))/(Time(i)Time(i1)); Velo_time(i1)=Time(i1)+0.5*(Time(i)Time(i1)); end >> figure >> plot(Velo_time, Velocity) >> title('Velocity vs. Time'),xlabel('Time'),ylabel('Velocity') Assignment 12 7 E7 Assignment 12 8 E7 University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz c) Use polyfit/polyder differentiation to compute acceleration and plot it.
>> hold on >> [out]=polyfit(Velo_time,Velocity,11); >> xxx=linspace(Velo_time(1),Velo_time(end),50); >> yyy=polyval(out,xxx) >> plot(xxx,yyy,'og') >> legend('velocity','fit to velocity') >> figure >> accel_time=linspace(Velo_time(1),Velo_time(end),50); >> [coeff]=polyder(out); >> acceleration=polyval(coeff,accel_time); >> plot(accel_time,acceleration) >> title('Acceleration vs. Time'),xlabel('Time'),ylabel('Acceleration') a) Write a function PGrowth that takes as its input arguments: P, the population,
r, the proportionality constant, and K, the carrying capacity, and computes the righthand side of Eq. (1). Solution: function Pdot = PGrowth(P,r,K) Pdot = r*P*(1  P/K); b) Write a function called PopulationGrowthODE45 that performs as follows: i) Takes as input arguments, P0 (the initial population), r (the proportionality constant), K (the carrying capacity), t0 (the initial time), and tf (the final time). (1) ii) Numerically solves the ODE in Eq. (1) from the initial time t0, to the final time tf, with the initial condition P0, using the matlab function 5. One equation for modeling population growth states that the change in population dP/dt can be expressed by the following ODE: dP P⎞ ⎛ = rP⎜1 − ⎟ dt K⎠ ⎝ This means that the rate of population change is dependent on the size of the existing population, P, and the amount of free resources available. In this equation, t is the time, P is the population at time t, r is the proportionality constant, and K is the maximum capacity of the system, often called the carrying capacity. This is equation is known as the Logistic Equation. ode45. The first input argument to ode45 should be the anonymous function handle PopulationGrowth = @(t,P) PGrowth(P,r,K) (See lecture 26 on Wednesday Nov. 26 for details). The solution of the ODE (1) is given by iii) K P0 e r t P(t ) = K + P0 (e r t − 1) where P0 = P(0) is the initial condition. (2) Plots the solution of the ODE in Eq. (1) (in blue) and the exact solution Eq. (2) (in green) as a function of time, and includes all the appropriate labels (e.g., title, axes, legends). Note that PopulationGrowthODE45 does not have any output arguments. Assignment 12 9 E7 Assignment 12 10 E7 University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz Run PopulationGrowthODE45 with the following input arguments: P0 = 10, r = 1.05, K = 2000, t0 = 0, tf = 15. At roughly what time does the system reach its carrying capacity? Do you believe that this is a realistic model for population growth? Why or why not? Solution: The function should return 2 output arguments, T and Y, just as ode45 does. The function should implement the Euler integration algorithm, as outlined in the Lecture 27 on Monday, December 1 lecture notes. Solution:
function [ T , Y ] = odeeuler( fh, tspan, y0, h) % uses same first 3 input arguments as ode45 function PopulationGrowthODE45(P0,r,K,t0,tf) PopulationGrowth = @(t,P) PGrowth(P,r,K); [T,P] = ode45(PopulationGrowth,[t0,tf],P0); Pexact=K*P0*exp(r*T)./(K+P0*(exp(r*T)1)); figure(1) plot(T,P,'b',T,Pexact,'r') % plus integration step size h n = length(y0); % assume yo is a column vector and determine % the order n (dimension of y0) T = (tspan(1):h:tspan(end))'; % define T (time) column vector M = length(T); % determine number of iterations use same convention as % ode45 to store results for Y >> PopulationGrowthODE45(10,1.05,2000,0,15) Y = zeros(M,n); Y(1,:) = y0'; Answer: This is not a realistic model, because it does not include death of any kind. It does not include interactions with predators. The carrying capacity K is a totally theoretical construct. No indication of change of seasons or breeding periods. Basically, this is a very contrived and theoretical approach to population models. 6. Write a function called odeeuler.
tspan, y0, h ) % perform Euler’s numerical integration algorithm for k=1:M1 Y(k+1,:) = Y(k,:) + h* fh(T(k), Y(k,:)')'; end function [ T , Y ] = odeeuler( fhan, 7. Write a function called PopulationGrowthEuler that The four input arguments are: i) • fhan: function handle that evaluates the righthand side of the ODE • tspan =[t0,tfinal]: simulation time span • y0: initial condition • h: fixed integration stepsize. ii) Takes as input arguments, P0 (the initial population), r (the proportionality constant), K (the carrying capacity), t0 (the initial time), the final time tf and the integration step size h. Numerically solves the ODE in Eq. (1), from the initial time t0, to the final time tf, with the initial condition P0, using your function Assignment 12 11 E7 Assignment 12 12 E7 University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz odeeuler. The first input argument to odeeuler should be the function [ T , Y ] = odemodeuler( fhan, tspan, y0, h) anonymous function handle
PopulationGrowth = @(t,P) PGrowth(P,r,K) that has the same input and output arguments as the function odeeuler. The function should implement the modified Euler algorithm, as outlined in the Lecture 27 on Monday, December 1 lecture notes. Solution: iii) Plots the solution of the ODE in Eq. (1) (in black) and the exact solution Eq. (2) (in green) as a function of time, and includes all the appropriate labels (e.g., title, axes, legends). function [ T , Y ] = odemodeuler( fh, tspan, y0, h) Run PopulationGrowthEuler three times each with the following input arguments
P0 = 10, r = 1.05, K = 2000, t0 = 0, tf = 15, and with the following % uses same first 3 input arguments as ode45 % plus integration step size h integration step sizes respectively: h = 1, .01, and .001. What do you notice about the difference between the exact solution and the Euler Method solution vs. the size of the integration step size, h? Solution: h2=h/2; n = length(y0); % assume % integration step size h yo is a column vector and determine % the order n (dimension of y0) T = (tspan(1):h:tspan(end))'; M = length(T); % define T (time) column vector % determine number of iterations % use same convention as ode45 to store results for Y, function PopulationGrowthEuler(P0,r,K,t0,tf,h) PopulationGrowth = @(t,P) PGrowth(P,r,K); [T,P] = odeeuler(PopulationGrowth,[t0,tf],P0,h); Y = zeros(M,n); Y(1,:) = y0'; % store results rowwise % perform Modified Euler’s numerical integration algorithm Pexact=K*P0*exp(r*T)./(K+P0*(exp(r*T)1)); figure(2) plot(T,P,'b',T,Pexact,'r') for k=1:M1 f1= fh(T(k), Y(k,:)); yp= Y(k,:)' + h*f1; % evaluate the first function call % do an Euler prediction >> for h = [1 0.01 0.001] figure PopulationGrowthEuler(10,1.05,2000,0,15,h) end % Do a secondorder approximation update Y(k+1,:) = Y(k,:) + h2* ( fh(T(k+1), yp)+ f1 )'; end Answer: The Euler solution converges to the exact solution for small h. However, it should be noted that this is not always the case. 8. Write a function called odemodeuler 9. Write a function called PopulationGrowthModEuler that takes the same input arguments and produces the same plots as PopulationGrowthEuler, but numerically solves the ODE in Eq. (1) using your funtion odemodeuler instead of
odeeuler. Assignment 12 13 E7 Assignment 12 14 E7 University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz Run PopulationGrowthModEuler three times each with the following input arguments P0 = 10, r = 1.05, K = 2000, t0 = 0, tf = 15, and with the following integration step sizes: h = 1, .01, and .001. Write a brief paragraph comparing the accuracies of the Euler and Modified Euler ODE integration algorithms for different step sizes. Solution: where t is the time, x, y, and z are the state variables, and σ, β, and ρ are scalar parameters. σ is called the Prandtl number and ρ is called the Rayleigh number. Usually σ = 10, β = 8/3 and ρ is varied. The system exhibits chaotic behavior for ρ = 28 but displays knotted periodic orbits for other values of ρ.
& The ODEs that describe the Lorentz attractor can be written in vector form Y = f (Y ) function = PopulationGrowthModEuler(P0,r,K,t0,tf,h) PopulationGrowth = @(t,P) PGrowth(P,r,K); [T,P] = odemodeuler(PopulationGrowth,[t0,tf],P0,h); Pexact=K*P0*exp(r*T)./(K+P0*(exp(r*T)1)); figure(3) plot(T,P,'b',T,Pexact,'r') function ⎡ Y1 ⎤ ⎡ σ ⋅ (Y2 − Y1 ) ⎤ d⎢ ⎥ ⎢ Y2 = Y1 ⋅ ( ρ − Y3 ) − Y2 ⎥ ⎥ dt ⎢ ⎥ ⎢ ⎢Y3 ⎥ ⎢ Y1 ⋅ Y2 − β ⋅ Y3 ⎥ ⎣⎦⎣ ⎦ (3) a) Write a function named LorenzODE dotY = LorenzODE(Y, sigma, beta, rho) >> for h = [1 0.01 0.001] figure PopulationGrowthModEuler(10,1.05,2000,0,15,h) end whose first input argument Y is a 3×1 column vector and sigma, beta, and rho are constants, and returns the right hand side of Eq. (3) (also a 3×1 column vector). Solution: Answer: The Modified Euler has superior performance characteristics than the Euler Method. ie. It gives more accurate results for the same step size. 10. The Lorenz attractor, named after Edward N. Lorenz, was derived as a simplified model of convection currents in the atmosphere. However, it is know for having chaotic, unpredictable behavior. It can be defined with the following three coupled differential equations: function dotY = LorenzODE(Y,sigma, beta,rho) dotY = zeros(3,1); dotY(1) = sigma*(Y(2)  Y(1)); dotY(2) = Y(1)*(rho  Y(3))  Y(2); dotY(3) = Y(1)*Y(2)  beta*Y(3); b) Use ode45 to numerically integrate the ODE Eq. (3) for: dy = x( ρ − z ) − y dt dz = xy − β z dt • • • σ = 10, β = 8/3, and ρ = 14
tspan = [0,20] dx = σ ( y − x) dt and the following set of initial conditions: Assignment 12 15 E7 Assignment 12 16 E7 University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz University of California, Berkeley Fall Semester 2008 Department of Mechanical Engineering Instructors: M. Frenklach, R. Horowitz i) ii) Y0 = [0;1;1]; ii) Y0 = [1;1;1]; iii) Y0 = [1;1;1]; Y0 = [1;1;1]; v) Y0 = [1;1;10]; c) Use ode45 to numerically integrate the ODE Eq. (3) for: • • • σ = 10, β = 8/3, and ρ = 28.
tspan = [0,20] Remember that first input argument to ode45 should be the anonymous function handle, such as and the following set of initial conditions:
i) Y0 = [0;1;1]; ii) Y0 = [0.001;1;1]; Lorenz = @(t,Y) LorenzODE(Y, sigma, beta, rho); Use the commands: Use the commands:
plot3(Y(:,1),Y(:,2),Y(:,3),’color’) , view(0,0), hold on plot3(Y(:,1),Y(:,2),Y(:,3),’color’) , view(0,0), hold on to do 3D plots of the solution (Y1 (t ), Y2 (t ), Y3 (t )) for the two initial conditions on the to do 3D plots the solution (Y1 (t ), Y2 (t ), Y3 (t )) for all five initial conditions on the same figure (using different values of color). Write a brief paragraph describing the results that you obtain.
>> sigma = 10; >> beta = 8/3; >> rho = 14; >> [T1,Y1] = ode45(@(t,Y) LorenzODE(Y,sigma,beta,rho),[0,20],[0;1;1]); >> [T2,Y2] = ode45(@(t,Y) LorenzODE(Y,sigma,beta,rho),[0,20],[1;1;1]); >> [T3,Y3] = ode45(@(t,Y) LorenzODE(Y,sigma,beta,rho),[0,20],[1;1;1]); >> [T4,Y4]=ode45(@(t,Y) LorenzODE(Y,sigma,beta,rho),[0,20],[1;1;1]); >> [T5,Y5]=ode45(@(t,Y) LorenzODE(Y,sigma,beta,rho),[0,20],[1;1;10]); same figure (using different values of color).
>> sigma = 10; >> beta = 8/3; >> rho = 14; >> [T6,Y6] = ode45(@(t,Y) LorenzODE(Y,sigma,beta,rho),[0,20],[0;1;1]); >> [T7,Y7] = ode45(@(t,Y) LorenzODE(Y,sigma,beta,rho),[0,20],... [0.001;1;1]); >> plot3(Y6(:,1),Y6(:,2),Y6(:,3),'r'),view(0,0),hold on >> plot3(Y7(:,1),Y7(:,2),Y7(:,3),'g') >> hold off d) Repeat all of part c) only that this time use tspan = [0,50]. Write a brief
>> plot3(Y1(:,1),Y1(:,2),Y1(:,3),'r'),view(0,0),hold on >> plot3(Y2(:,1),Y2(:,2),Y2(:,3),'g') >> plot3(Y3(:,1),Y3(:,2),Y3(:,3),'b') >> plot3(Y4(:,1),Y4(:,2),Y4(:,3),'c') >> plot3(Y5(:,1),Y5(:,2),Y5(:,3),'y') >> hold off paragraph describing the results that you obtain. Answer: Even for a longer time span, the Lorenz attractor for the two initial starting points (that happen to be very close together) stay together. Assignment 12 17 E7 Assignment 12 18 E7 ...
View
Full
Document
This note was uploaded on 11/01/2009 for the course ENGLISH 7 taught by Professor Sengupta during the Spring '09 term at Berkeley.
 Spring '09
 Sengupta

Click to edit the document details