Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
NUMERICAL_METHODS_COMPILATION_IN_MATLAB
Typology: Thesis
1 / 75
This page cannot be seen from the preview
Don't miss anything!
PHYSICISTS
ENRIQUE ALFONSO FONSECA MANRIQUE
Universidad de Guadalajara Centro de Ciencias Exactas e Ingenierías
P L 2017 GRAND CANON
2
The plot command is by far the most simple way to create a two-dimensional (2- D) plot of the data in y versus the corresponding values of x. The general syntax is as follows
plot(x,y, ’LineSpec’, ’PropertyName’, PropertyValue)
where x and y can be either scalars or vectors^1 ; ,“LineSpec” (which is typed as strings) formats your line (color, type of line and markers) and“PropertyName” and ”PropertyValue” are properties with values that define, for example, the line width, marker’s size, and other geometrical attributes. We will review briefly both of these options.
At this point, we should have an intuitive notion of how to use the plot command and hereunder we shall use some of the aformentioned properties and introduce important details about plotting.
clear ; clc
x = linspace(0,2pi,200); y = sin(x).ˆ2.cos(2x); yd = 2sin(x).cos(x).cos(2x) - 2sin(x).ˆ2.sin(2x);
figure(’Units’,’inches’,... ’Position’,[3 3 8 5],... ’PaperpositionMode’,’auto’);
plot(x, y,’-r’, x, yd, ’--b’)
set(gca,... ’Units’,’normalized’,... ’YTick’, -2:.5:2,... ’XTick’, -2:1:10,... ’Position’,[.15 .2 .75 .7],... ’FontUnits’,’points’,... ’FontWeight’,’normal’,... ’FontSize’,11,... ’FontName’,’Times’)
grid on
(^1) It is important to mention that if both x and y are vectors or matrices they must be equal in size; if one of them is a scalar and the other either a scalar or vector then the plot function creates a visual array of discrete points which must be visualized via a marker (for later reference see ??)
Evaluation of Functions and Graphing 3
xlabel(’x’) ylabel(’y’) title(’Nice Plot’) legend({’f(x)’,’First Derivative’},’FontSize’, 12,’Location’,’southeast’)
print(’nice’,’-depsc’)
0 1 2 3 4 5 6 7 x
-1.
-0.
0
1
2
y
Nice Plot
f(x) First Derivative
Fig. 1: A nice graph created by special settings in our script. Without it, our plot would have a poor resolution.
We have used some extra commands which haven’t been covered previously. First of all, when we export a MatLab figure and we want to use it on a document we’ll find some nasty results; mainly, the figure won’t have a nice resolution and won’t be informative, which is a serious issue. For this reason, we have to print our figure directly from the command windows or in our script. First, the figure command sets the size of the image. Next, gca configures the ticks our axis. Finally, the print command saves our figure into a desired format. The result of this procedure is illustrated in Fig.1.1.
Evaluation of Functions and Graphing 5
Fig. 2: Different types of settings
where the density of the grid is determined by the number of elements in the vectors. For example,
x = -1:3; y = 1:4; [X,Y] = meshgrid(x,y)
X =
-1 0 1 2 3 -1 0 1 2 3 -1 0 1 2 3 -1 0 1 2 3
Y =
1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4
where we can clearly see the coordinates.
6
Graphing surfaces in 3-D space by hand sometimes may be a little too diffcult be- cause the functions which we want to graph are really complex. Nevertheless, there are some alternative tools which can help us visualize the graph of a function f (x, y). The level curves of a function f (x, y) are the curves of the equation with the form f (x, y) = k where k is a constant. A set of level curves z = f (x, y) is called a contour map. This method provides us with a topographic projection of the equation we are studying. The basic syntax is:
contour(X,Y,Z,n) contour(X,Y,Z,v) contour(...,LineSpec)
where the contour plot displays isolines of matrix Z, n indicates the number of contour levels and v is a monotonically increasing vector 2.
In this section we will introduce two important tools in programming known as con- ditionals and loops. We’ll also apply some of the methods and commands that were mentioned in the previous section aiming to illustrate all of them with an introduc- tory problem in classical mechanics: projectile motion.
A relational operator compares two numbers and determines if the built-in requisite is fulfilled, this is either true (1) or false (0). Logical operators are concerned with statements which might be true or false and, correspondingly, assigns a (1) or (0). These operators (in combinations with other commands) are used to “orquestrate” the flow of the computer program.
Relational Operators
Relational operators are part of a programming language that is used to test or define a certain condition; most of these conditions are numerical equalities or inequalities, such as, 5 = 5 or 6 > 3. These operators vary from language to language; those used in MatLab are described in Table 1.1.
(^2) This is particularly useful when we desire to display a countour level at a particular value.
8
In programming, loops specify iterations, which allows code to be executed repeat- edly. MatLab has two basic loops: for loop and the while loop. The first of these has the basic syntax:
for i = first to last -- statement -- end
where, it is necessary to know how many iterations will be used. On the other hand, a while loop has the basic syntax:
while (boolean condition) then: -- statement -- end
combining a condition and an iteration, allowing the code to be executed repeatedly based on a Boolean condition.
In this section we’ll use some of the previews methods and tools to solve graphically and numerically a simple equation known as ”projectile motion”. First, let’s explain some of the physics involved during this 2-D motion.
Basic Equations
Consider we’re analyzing the motion of a baseball under the influence of a gravita- tional field. To describe this motion we define the vector position r(t) and velocity vector as:
dv dt
= −g yˆ ,
dr dt
= v,
where we’re not considering wind resistance. Notice that in the previous equation the acceleration during the motion is only observerd in the ˆy direction. This means we could express the coordinate y as a function of time:
y(t) = yo + vosin(θ )t −
gt^2 (1)
where we have denoted yo (t=0) as the initial vertical position, vo (t=0) initial ve- locity and θ is the angle above the horizontal at which velocity is described. It is important to mention that these equations only hold when there is no air resistance.
Evaluation of Functions and Graphing 9
The velocity in the y-component is described as vosinθ , which we shall refer to it as voy from hereon. Geometrically, Eq. 1,1 describes a parabola that starts when y 0 at t = 0 and should end at y = 0 at some arbitrary time t 3. This is a physical constraint because it doesn’t allow the particle to go down Earth into oblivion, such a case is most un- desirable. However, we could consider that the particle is bouncing, which may be a reasonable case for our baseball ball, nevertheless, we will analyze this fact later on. It is important to know how a function is changing through space, for this specific motion we’re interested in finding a maximum, that is the maximum height reached. Physically, this happens when the velocity is equal to 0 and it represents the point of space in which velocity changes sign, from positive to negative and starts falling. It is posible to find this maximum using the derivative of the function,
y′(t) = voy − gt, tmax =
voy g
determines the time at which the maximum ocurrs, and we shall denote the maxi- mum height as ymax.
Simple Projectile Motion [3][4]
The procedure is quite straight forward: to create a plot of the trajectory of the baseball ball we should choose an arbitrary length for the time t, define the initial conditions yo, voy and write down the equation. For example:
clear ; clc
t = 0:.01:4; %Define interval of time
y0 = 0 ; %Initial position v0 = 20; %Initial velocity g = -9.8;
%Trajectory equation
y = y0 + v0t + 0.5g*t.ˆ2;
%Plot height versus time elapsed
plot(t,y) grid on
This is obviously a very primitive plot which provides little information about what is going on during the process of falling. We could enhance this plot by including the
(^3) The time at which y = 0 (when yo 6 = 0) is the positive root of the quadratic equation y = −voy±√v^2 oy+ 2 yog −g , evidently because time can’t be negative.
Evaluation of Functions and Graphing 11
1 until it reaches 3. The velocity in function of k is rather derived arbitrarily, thus, there exists more ways in that it could be established. Once an initial velocity is determined the program continues to use Eqs. (1.1) and (1.2) to calculate the trajec- tory and velocity correspondingly. Then, it proceeds to graph both functions in the same window, but in separate plots which is accomplished by using the command subplot. At this point the program restarts the loop but now uses the value of k=2, proceeds as described previously, stores the data and plots it on the same graphs.
12
Bouncing Ball
Previously, we weren’t considering the constraints of the problem. A constraint is a geometrical limitation, in projectile motion is quite obvious that these limitation is when the ball hits the floor which causes the ball to bounce and deliver some energy to the ground provoking a shorter maximum height until it stays on the ground. Nevertheless, we’ll focus on the ideal case where there is no energy lost by the ball so the energy is conserved along the process thus, the maximum height should stay the same. In this example we used another way to create an interval using a for loop. For our time interval we defined a step size which determines how fast it will change through iterations. Thus, for some k it will correspond some value of t being time a function of k. It is necesary to redifine our initial conditions for every value of k. The second important moment is defining the instant at which the ball hits the floor. At such instant the ball will change its velocity from negative to positive and will now travel into air again. Then, the constraint is that y ≥ 0. This, translated into MatLab language, can ba stated as: if y < 0 the ball should change from v to −v (being −v > 0) as seen in Figure 1.4 in which the velocity has a sudden change from negative to positive. It
Fig. 4: It is observed a quite obvious discontinuity in the velocity because of the sudden change of velocity upon collision.
14
Code for ”Several Initial Conditions”
clear ; clc
%Projectile motion with different initial conditions
t = (0:.1:5); g = -9.81; y0 = 20; %(m)
% Construction of different velocities using a "for" loop for k = 1:
v0 = 20*(k-2); %(m/s)
y (k, :) = 0.5gt.ˆ2 +v0.t +y0; %Index k yd(k, :) = gt+ v0;
end
subplot(2,1,1) plot(t,y, ’--’, ’LineWidth’, 1.5) hold on plot([0, 5],[0 0],’k’) grid on
title(’Trajectory’) xlabel(’Time (s)’) ylabel(’Position (m)’) legend(’v_o = -20m/s’,’v_o = 0m/s’,’v_o = 20m/s’) axis([0 5 -30 50])
subplot(2,1,2) plot(t,yd) grid on
title(’Velocity’) xlabel(’Time(s)’) ylabel(’Velocity (m/s)’) legend(’v_o = -20m/s’,’v_o = 0m/s’,’v_o = 20m/s’)
Evaluation of Functions and Graphing 15
Code for ”Boucing Ball”
clear ; clc
%Physical Parameters g = -9.8;
%Step size dt = 0.01;
%Initial Conditions
yo=20; % (m) vo=20; % (msˆ-1) clf
for k=1: t = dt(k-1); y = 0.5gdt.ˆ2+vodt+yo; v = g*dt+vo;
yo=y; vo=v;
if y< vo=-v; y=0; end
subplot(2,1,1); plot(t,y,’ok’,’MarkerSize’,0.5) hold on plot([0 14], [0 0],’k’); xlabel(’Time (s)’) ylabel(’Height (m)’) grid on axis([0 14 -10 50]) title(’Trajectory’)
subplot(2,1,2); plot(t,v,’ok’,’MarkerSize’,0.5)
xlabel(’Time (s)’) ylabel(’Velocity (m/s)’) grid on axis([0 14 -30 20]) title(’Velocity (m/s)’) hold on
end
2
Fig. 1 The interval where the real root lies (XNS) has a change of sign and there- fore satifies the condition f (a) f (b) < 0. Bisection method is only defined for functions which are continous in the interval [a, b].
Mathematically, we search for a value of x such that the function f (x) = 0 in a given interval [a, b] satisfying f (a) f (b) < 0 implying a change of sign in the function lying on the interval. Therefore, the first estimate of the solution is
xs 1 = a + b 2
Now, it is necessary to determine whether the true solution lies between [a, xs 1 ] or [xs 1 , b], to check these conditions we use
f (a) f (xs 1 ) < 0 (3) f (a) f (xs 1 ) > 0 (4)
if Eq(1.3) is true, the solution lies in the interval [a, xs 1 ]; otherwise, use the interval [xs 1 , b]. Once selected, return to Eq(1.2) and replace the new limit and evaluate a new solution. We define the stop point of the iterations once the following condition is satisfied
Tolerance < |b − a| 2
where the Tolerance is defined by the user. Let us now show an example applying the bisection method.
The code described in the Appendix shows a new tool in MatLab called anonymous function which is a very practical choice when a relatively simple mathematical expression is used many times in a single code. In our case, we define an arbitrary non-linear function. We chose
F = @(x) x.ˆ3 -3.x.ˆ2+3x - sin(x) + 2;
Finding Real Zeros 3
as the function we want to solve. Notice that the function is named F and the argu- ment (independent variables) is enclosed by @(), in general, an anonymous function is created as:
Name = @ (x1,x2,...,xn) expression
where expression refers to the mathematical expression.
-1 -0.8 -0.6 -0.4 -0.2 0 x
0
1
2
3
4 Function a b F(c)
Our first interval [a, b] is fixed and we specify a tolerance tol as well as a maximum of iterations max in the for loop. Before the for loop we wish that our interval satifies that F(a)F(b) < 0. If the condition is not met an error is displayed which warns the user of this fact. We added, for the sake of utility, a table in which various variables are displayed for each value of k. The rest of the code shows the procedure described in Section 1.1.
The secant method[2] is a scheme for finding a numerical solution of an equation of the form f (x) = 0. It consists of using two points (x 1 , x 2 ) in a neighborhood of the solution [x 1 , x 2 ] to define secant line and the point of intersection x 3 is the new estimate of the solution. Using Figure 1.2b, the slope of the secant line is given by
