Matlab Workbook, Exams of Differential Equations

Download Matlab Workbook and more Exams Differential Equations in PDF only on Docsity!

MATLAB WORKBOOK

CME 102 Winter 2008-

Eric Darve

Hung Le

2/55 CME 102 Matlab Workbook 2008-

Introduction

This workbook aims to teach you Matlab and facilitate the successful integration of Matlab into the CME 102 (Ordinary Differential Equations for Engineers) curriculum. The workbook comprises three main divisions; Matlab Basics, Matlab Programming and Numerical Methods for Solving ODEs. These divisions are further subdivided into sections, that cover specific topics in Matlab. Each section begins with a listing of Matlab commands involved (printed in boldface), continues with an example that illustrates how those commands are used, and ends with practice problems for you to solve. The following are a few guidelines to keep in mind as you work through the examples:

a) You must turn in all Matlab code that you write to solve the given problems. A convenient method is to copy and paste the code into a word processor.

b) When generating plots, make sure to create titles and to label the axes. Also, include a legend if multiple curves appear on the same plot.

c) Comment on Matlab code that exceeds a few lines in length. For instance, if you are defining an ODE using a Matlab function,explain the inputs and outputs of the function. Also, include in-line comments to clarify complicated lines of code.

Good luck!

4/55 CME 102 Matlab Workbook 2008-

Solution:

a) >> A = zeros(2,4) A = 0 0 0 0 0 0 0 0

b) >> L = 1 : 2 : 21 L = 1 3 5 7 9 11 13 15 17 19 21

c) >> S = sum(L) S = 121

d) >> A = [2, 3, 2; 1 0 1] A = 2 3 2 1 0 1

1.1.2 Your Turn

a) Create a 6 x 1 vector a of zeros using the zeros command.

b) Create a row vector b from 325 to 405 with an interval of 20.

c) Use sum to find the sum a of vector b’s elements.

1.2 Matrix and Vector Operations

Commands:

• Element-by-element addition. (Dimensions must agree)

• Element-by-element subtraction. (Dimensions must agree) .* Element-by-element multiplication. (Dimensions must agree) ./ Element-by-element division. (Dimensions must agree) .^ Element-by-element exponentiation. : When used as the index of a matrix, denotes “ALL” elements of that dimension. A(:,j) j-th column of matrix A (column vector). A(i,:) i-th row of matrix A (row vector). .’ Transpose (Reverses columns and rows). ’ Conjugate transpose (Reverses columns and rows and takes complex conjugates of elements).

• Matrix multiplication, Cayley product (row-by-column, not element- by-element).

CME 102 Matlab Workbook 2008-2009 5/

1.2.1 Example

a) Create two different vectors of the same length and add them.

b) Now subtract them.

c) Perform element-by-element multiplication on them.

d) Perform element-by-element division on them.

e) Raise one of the vectors to the second power.

f) Create a 3 × 3 matrix and display the first row of and the second column on the screen.

Solution:

a) >> a = [2, 1, 3]; b = [4 2 1]; c = a + b

c =

6 3 4

b) >> c = a - b

c =

-2 -1 2

c) >> c = a .* b

c =

8 2 3

d) >> c = a ./ b

c =

0.5000 0.5000 3.

e) >> c = a .^ 2

c =

4 1 9

f) >> d = [1 2 3; 2 3 4; 4 5 6]; d(1,:), d(:,2)

ans =

1 2 3

CME 102 Matlab Workbook 2008-2009 7/

ylabel(’...’) Similar to the xlabel command. grid on Adds grid lines to the current axes. grid off Removes grid lines from the current axes.

1.3.1 Example

Let us plot projectile trajectories using equations for ideal projectile motion:

y(t) = y 0 −

gt^2 + (v 0 sin(θ 0 )) t, x(t) = x 0 + (v 0 cos(θ 0 )) t,

where y(t) is the vertical distance and x(t) is the horizontal distance traveled by the projectile in metres, g is the acceleration due to Earth’s gravity = 9.8 m/s^2 and t is time in seconds. Let us assume that the initial velocity of the projectile v 0 = 50.75 m/s and the projectile’s launching angle θ 0 = 512 π radians. The initial vertical and horizontal positions of the projectile are given by y 0 = 0 m and x 0 = 0 m. Let us now plot y vs. t and x vs. t in two separate graphs with the vector: t=0:0.1:10 representing time in seconds. Give appropriate titles to the graphs and label the axes. Make sure the grid lines are visible.

Solution:

We first plot x and y in separate figures:

t = 0 : 0.1 : 10; g = 9.8; v0 = 50.75; theta0 = 5pi/12; y0 = 0; x0 = 0; y = y0 - 0.5 * g * t.^2 + v0sin(theta0).t; x = x0 + v0cos(theta0).*t;

figure; plot(t,x); title(’x(t) vs. t’); xlabel(’Time (s)’); ylabel(’Horizontal Distance (m)’); grid on;

figure; plot(t,y); title(’y(t) vs. t’); xlabel(’Time (s)’); ylabel(’Altitude (m)’); grid on;

8/55 CME 102 Matlab Workbook 2008-

0 5 10 0

50

100

150

x(t) vs. t

Time (s)

Horizontal Distance (m) 0 5 10

0

50

100

150

y(t) vs. t

Time (s)

Altitude (m)

1.3.2 Your Turn

The range of the projectile is the distance from the origin to the point of impact on horizontal ground. It is given by R = v 0 cos(θ 0 ). To estimate the range, your trajectory plots should be altered to have the horizontal distance on the x-axis and the altitude on the y-axis. This representation will clearly show the path of the projectile launched with a certain initial angle. This means you will have to plot y vs. x. Observing the formula for the projectile’s range, we see that to increase the range we will have to adjust the launching angle. Use the following adjusted angles to create two more trajectory plots (y vs. x), one for each angle, and determine which launching angle results in a greater range:

θ^10 =

5 π 12

radians and

θ^20 =

5 π 12

radians.

The time vectors for these angles should be defined as t = 0:0.1:9 and t = 0:0.1:8 respectively.

1.4 Plotting Multiple Functions I

Commands:

plot(x,y) Creates a plot of y vs. x. plot(x,y1,x,y2, ...)

Creates a multiple plot of y1 vs. x, y2 vs. x and so on, on the same fig- ure. MATLAB cycles through a predefined set of colors to distinguish between the multiple plots. hold on This is used to add plots to an existing graph. When hold is set to on, MATLAB does not reset the current figure and any further plots are drawn in the current figure. hold off This stops plotting on the same figure and resets axes properties to their default values before drawing a new plot. legend Adds a legend to an existing figure to distinguish between the plotted curves. ezplot(’f(x)’, [x0,xn])

Plots the function represented by the string f(x) in the interval x0 ≤ x ≤ xn.

10/55 CME 102 Matlab Workbook 2008-

0 1 2 3 4 5 6

−0.

0

x

High Frequency Cosine Function

y

1.4.2 Your Turn

a) Using the plot command for multiple plots, plot y = atan(x) and y = acot(x) on the same graph for values of x defined by x = -pi/2:pi/30:pi/2.

b) Using the plot command for a single plot and the hold commands, plot y = atan(x) and y = acot(x) on the same graph for values of x defined by x = -pi/2:pi/30:pi/2.

c) Using the ezplot command, plot y =^23 sin(9πx), for values of x such that 0 ≤ x ≤ 2 ∗ pi.

1.5 Plotting Functions II

Commands:

log(n) Calculates the natural logarithm (base e) of n. semilogy(x,y) Graphs a plot of y vs. x using a logarithmic scale (powers of ten) on the y-axis. semilogx(x,y) Graphs a plot of y vs. x using a logarithmic scale (powers of ten) on the x-axis. loglog(x,y) Graphs a plot of y vs. x using a logarithmic scale (powers of ten) on both axes. The logarithmic scales prove most useful when the value spans multiple orders of magnitude.

1.5.1 Example

Graph the efficiency of several programming algorithms according to big-O notation, a method of describing the running time of algorithms. Each expression represents the scale by which an algorithm’s computation time increases as the number of its input elements increases. For example, O(n) represents an algorithm that scales linearly, so that its computation time increases at the same rate as the number of elements. The algorithms you must graph have the following big-O

CME 102 Matlab Workbook 2008-2009 11/

characteristics:

Algorithm #1: O(n) Algorithm #2: O(n^2 ) Algorithm #3: O(n^3 ) Algorithm #4: O(2n) Algorithm #5: O(en)

After generating an initial graph with ranging from 0 to 8, use logarithmic scaling on the y-axis of a second graph to make it more readable. You can also use the mouse to change the y-axis scale. Go to the main menu of the figure, click Edit>Axes Properties... , the property editor dia- logue will pop out. There, you can also change the font, the range of the axes,... Try to play with it.

Solution:

n=0:0.01:8; plot(n,n,n,n.^2,n,n.^3,n,2.^n,n,exp(n)) title(’Big-O characteristics of Algorithms: Linear Plot’) ylabel(’Estimate of Running Time’) xlabel(’n (number of elements)’) legend(’O(n)’,’O(n^2)’,’O(n^3)’, ’O (2^n)’,’O(e^n)’) grid on;

(^00 2 4 6 )

500

1000

1500

2000

2500

3000

Big−O characteristics of Algorithms: Linear Plot

Estimate of Running Time

n (number of elements)

O(n) O(n^2 ) O(n^3 ) O (2n) O(en)

n = 0:0.01:8; semilogy(n,n,’b’,n,n.^2,’r’,n,n.^3,’g’,n,2.^n,’c’,n,exp(n),’k’) title(’Big-O characteristics: Logarithmic Plot’) ylabel(’Estimate of Running Time’) xlabel(’n (number of elements)’) legend(’O(n)’,’O(n^2)’,’O(n^3)’, ’O(2^n)’,’O(e^n)’)

CME 102 Matlab Workbook 2008-2009 13/

... The ellipses can be used to break up long lines by providing a contin- uation to the next line. Strings must be ended before the ellipses but can be immediately restarted on the next line. Examples below show this.

Note: Neglecting the command clear can cause errors because of previously defined variables in the workspace.

fprintf: This is an example of how to use fprintf to display text to the command window.

fprintf (’\nOrdinary Differential Equations are not so ordinary.\n’); fprintf (’-------------------------------------------------------’... ’----------------\n’); fprintf (’This course is CME %g: ODEs for Engineers. My expected’... ’ grade is %g\n’,102,100); x = 100; y = 96; fprintf (’The previous course was CME %g: Vector Calculus for ’... ’Engineers. My grade was: %g\n’,x,y);

The Matlab command window displays:

Ordinary Differential Equations are not so ordinary.

This course is CME 102: ODEs for Engineers. My expected grade is 100 The previous course was CME 100: Vector Calculus. My grade was: 96

The command fprintf takes a string and prints it as is. The character \n is one of several “Escape Characters” for fprintf that can be placed within strings given to fprintf. \n specifies a new line. %g is one of many “Specifiers” that fprintf uses and it represents a placeholder for a value given later in the call to fprintf. The order of the arguments given to fprintf determine which %g is replaced with which variable or number. Experiment with the code above to see what \n can do and how %g can be used.

M-Files/Scripts: Everything we have done so far has been in MATLABs interactive mode. However, MATLAB can execute commands stored in a regular text file. These files are called scripts or ’M-files’. Instead of writing the commands at the prompt, we write them in a script file and then simply type the name of the file at the prompt to execute the commands. It is almost always a good idea to work from scripts and modify them as you go instead of repeatedly typing everything at the command prompt. A new M-file can be created by clicking on the “New M-file” icon on the top left of the Command Window. An M-file has a .m extension. The name of the file should not conflict with any existing MATLAB command or variable. Note that to execute a script in an M-file you must be in the directory containing that file. The

14/55 CME 102 Matlab Workbook 2008-

current directory is shown above the Command Window in the drop down menu. You can click on the “... ” icon, called “Browse for folder”, (on the right of the drop-down menu) to change the current directory. The % symbol tells MATLAB that the rest of the line is a comment. It is a good idea to use comments so you can remember what you did when you have to reuse an M-file (as will often happen). It is important to note that the script is executed in the same workspace memory as everything we do at the prompt. We are simply executing the commands from the script file as if we were typing them in the Command Window. The variables already existing before executing the script can be used by that script. Similarly, the variables in the script are available at the prompt after executing the script.

2.1.1 Example

After your 30 years of dedicated service as CEO, TI has transfered you to a subdivision in the Himalayas. Your task as head of the subdivision is to implement transcendental functions on the Himalayan computers. You decide to start with a trigonometric function, so you find the following Taylor Series approximation to represent one of these functions:

sin(x) = x − x^3 3!

x^5 5!

− · · · + (−1)n^ x^2 n+ (2n + 1)!

However, since the computers in the Himalayas are extremely slow (possibly due to the high alti- tudes), you must use the Taylor Series as efficiently as possible. In other words, you must use the smallest possible number of terms necessary to be within your allowed error, which is 0.001. You will use x = 3 as the value at which to evaluate the function.

a) Compute and display the exact error of using the first 2 and 5 terms of the series as compared to the actual solution when the function is evaluated at x = π 3.

b) Compute and display the number of terms necessary for the function to be within the allowed error.

Solution:

a) CODE from M-file

clear; x = pi/3;

% Iterate 2 terms for our approximation. SIN_Approx2 = 0; for j=0: SIN_Approx2 = SIN_Approx2 + (-1)^jx^(2j+1)/factorial(2*j+1); end SIN_Error2 = abs(SIN_Approx2 - sin(x));

% Iterate 5 terms for our approximation. SIN_Approx5 = 0; for j=0: SIN_Approx5 = SIN_Approx5 + (-1)^jx^(2j+1)/factorial(2*j+1);

16/55 CME 102 Matlab Workbook 2008-

2.1.2 Your Turn

The moment you finish implementing the trigonometric functions, you realize that you have for- gotten your favorite function: the exponential! The Taylor Series for the exponential function is

ex^ = 1 + x + x^2 2!

x^3 3!

xn n!

a) Using a for loop, compute and display the error of the Taylor series when using 12 and 15 terms, each as compared to the exponential function when evaluated at x = 2.

b) Using a while loop, compute and display the number of terms necessary for this function to have an error within your allowed range of 0.001 when evaluated at x = 3. Display this error.

2.2 Symbols and Strings: Differentiation and Integration

Commands:

syms x y z Declares variables x, y, and z as symbolic variables. diff(f) Symbolically computes derivative of the single-variable symbolic ex- pression f. diff(f,n) Symbolically computes the n-th derivative of the single-variable sym- bolic expression f. diff(f,’x’) Symbolically differentiates f with respect to x. diff(f,’x’,n) Symbolically computes the n-th derivative with respect to x of f. int(f) Symbolically integrates the single-variable symbolic expression f. int(f,x) Symbolically integrates f with respect to x.

2.2.1 Example

Matlab is capable of symbolically computing the derivatives and integrals of functions. This can be very convenient to check your pencil and paper solutions; in some cases Matlab may be able to integrate functions which require very difficult change of variables. We will also learn in the next section how Matlab can be used to symbolically solve differential equations. For now we focus on derivatives and integrals.

The first step is to define symbolic variables using the syms command. Note that symbolic variables are different from numerical variables. They can not assume any numerical value. They are simply symbols or letters that can be manipulated symbolically.

Once these variables are declared we can define symbolic functions, differentiate and integrate them as is shown in the following examples:

a) Define the following function using syms:

f (x) = x^2 ex^ − 5 x^3.

b) Compute the integral, and the first and second derivatives of the above function symbolically.

CME 102 Matlab Workbook 2008-2009 17/

c) Define the following function using syms:

f (x, y, z) = x^2 ey^ − 5 z^2.

Compute the integral with respect to x and second derivative with respect to z.

Solution:

a) >> syms x

f = x^2exp(x) - 5x^3;

b) >> int(f) ans = x^2exp(x)-2xexp(x)+2exp(x)-5/4*x^

diff(f) ans = 2xexp(x)+x^2exp(x)-15x^

diff(f,2) ans = 2exp(x)+4xexp(x)+x^2exp(x)-30*x

c) >> syms x y z

f = x^2exp(y) - 5z^2;

int(f,’x’) ans = 1/3x^3exp(y)-5z^2x

diff(f,’z’,2) ans =

2.2.2 Your Turn

a) Compute

(x − 1)(x − 3)

dx.

b) Compute

x^2 e−^3 xdx.

c) Compute the first and the second derivatives of atan

( (^) x 5

CME 102 Matlab Workbook 2008-2009 19/

c) Find the particular solution of the following ODE:

dy dt = e−y y(0) = 0

d) Find the particular solution of the following ODE:

dy dt

• Ky = A + B cos

πt

y(0) = 0

where A, B and K are constants.

Solution:

a) >> dsolve(’Dy+ytan(t) = sin(2t)’) ans = -2cos(t)^2+cos(t)C

b) >> dsolve(’D2y-Dy=1+tsin(t)’) ans = -1/2cos(t)-1/2tsin(t)-sin(t)+1/2cos(t)t+exp(t)*C1-t+C

c) >> dsolve(’Dy=exp(-y)’,’y(0)=0’) ans = log(t+1)

d) dsolve(’Dy+Ky = A+Bcos((1/12)pit)’,’y(0)=0’) ans =

exp(-Kt)(-144BK^2-Api^2-144AK^2)/(144K^3+Kpi^2) +(144AK^2+Api^2+144BK^2cos(1/12pit) +12Bpisin(1/12pit)K)/K/(144K^2+pi^2)

2.3.2 Your Turn

a) Find the general solution of the following first order ODE:

dy dt

− y = e^2 t.

b) Find the general solution of the following second order ODE:

d^2 y dt^2

dy dt

• 16y = 0.

c) Find the particular solution in a) given the following initial conditions:

y(0) = 5.

20/55 CME 102 Matlab Workbook 2008-

3 Forward Euler

Commands:

norm(v) Calculates the Euclidean norm of a vector v. length(v) Returns the length of vector v. axis equal Sets the aspect ratio so that each axis has the same scale.

Forward Euler presents the most basic way to solve a differential equation numerically. Even though it has the most stringent stability criterion and the lowest accuracy of all numerical solvers, it is the simplest and easiest to program. In this section, we will explore:

• The algorithm and code used to program Forward Euler
• Application of Forward Euler to solving differential equations
• Computing local and global error
• Forward Euler’s order of accuracy
• Forward Euler’s stability limitations

Algorithm for Forward Euler

yn+1 = yn + h · y′ n = yn + h · f (tn, yn)

The forward Euler algorithm is an explicit means of obtaining the value of y at each time step using the previous time step and the local slope or derivative. The equation is derived from approximating the derivative y′^ with the difference quotient

yn+1 − yn h ≈ y n′ = f (tn, yn).

As the boxed algorithm expresses, we solve for each subsequent time step yn+1 by adding the approximate change in y to the previous value of y. The last statement in the equation follows our convention of designating the slope or derivative as the function f (tn, yn). Thus, by decreasing the time step h by which we march across our interval, we can make our solution more and more accurate. This algorithm is explicit because all values for the next time step are readily isolable from previous known values, making evaluation simple and systematic.

Code for Forward Euler

To solve any ordinary differential equation that is expressible in the explicit form

y′ n = f (tn, yn),

we pass a Matlab function yprime.m that contains the code for computing f (tn, yn) to another function, say, forwardEuler.m. (The latter performs the numerical integration according to the forward Euler method.) Once we have the forward Euler algorithm implemented in Matlab, we