**UFBA**

# Partial Differential Equations in MATLAB

(Parte **1** de 4)

MATLAB Tutorial to accompany

Partial Differential Equations: Analytical and Numerical

Methods, 2nd edition by

Mark S. Gockenbach (SIAM, 2010)

MATLAB Tutorial | 1 |

Introduction | 3 |

About this tutorial | 3 |

About MATLAB | 3 |

MATLAB M-Book | 3 |

Getting help with MATLAB commands | 4 |

Getting started with MATLAB | 4 |

Vectors and matrices in MATLAB | 8 |

More about M-Books | 9 |

Simple graphics in MATLAB | 9 |

Symbolic computation in MATLAB | 15 |

Manipulating functions in MATLAB | 18 |

Saving your MATLAB session | 20 |

About the rest of this tutorial | 20 |

Chapter 1: Classification of differential equations | 21 |

Chapter 2: Models in one dimension | 24 |

Section 2.1: Heat flow in a bar; Fourier's Law | 24 |

Solving simple boundary value problems by integration | 25 |

The MATLAB solve command | 26 |

Chapter 3: Essential linear algebra | 31 |

Section 3.1 Linear systems as linear operator equations | 31 |

Section 3.2: Existence and uniqueness of solutions to Ax=b | 32 |

Section 3.3: Basis and dimension | 35 |

Symbolic linear algebra | 37 |

Programming in MATLAB, part I | 38 |

Defining a MATLAB function in an M-file | 38 |

Optional inputs with default values | 43 |

Comments in M-files | 44 |

M-files as scripts | 44 |

Section 3.4: Orthogonal bases and projection | 46 |

Working with the L2 inner product | 46 |

Section 3.5: Eigenvalues and eigenvectors of a symmetric matrix | 53 |

Numerical eigenvalues and eigenvectors | 53 |

Symbolic eigenvalues and eigenvectors | 53 |

Review: Functions in MATLAB | 56 |

Chapter 4: Essential ordinary differential equations | 59 |

Section 4.2: Solutions to some simple ODEs | 59 |

Second-order linear homogeneous ODEs with constant coefficients | 60 |

A special inhomogeneous second-order linear ODE | 61 |

First-order linear ODEs | 62 |

Section 4.3: Linear systems with constant coefficients | 64 |

Inhomogeneous systems and variation of parameters | 66 |

Programming in MATLAB, Part I | 68 |

Conditional execution | 69 |

Passing one function into another.............................................................................................. ........... 70

Section 4.4: Numerical methods for initial value problems | 72 |

Programming in MATLAB, part I | 78 |

Efficient MATLAB programming | 81 |

More about graphics in MATLAB | 81 |

Chapter 5: Boundary value problems in statics | 83 |

Section 5.2: Introduction to the spectral method; eigenfunctions | 83 |

Section 5.5: The Galerkin method | 88 |

Computing the stiffness matrix and load vector in loops | 92 |

Section 5.6: Piecewise polynomials and the finite element method | 93 |

Computing the load vector | 95 |

Computing the stiffness matrix | 96 |

The nonconstant coefficient case | 99 |

Creating a piecewise linear function from the nodal values | 100 |

More about sparse matrices | 101 |

Chapter 6: Heat flow and diffusion | 104 |

Section 6.1: Fourier series methods for the heat equation | 104 |

Section 6.4: Finite element methods for the heat equation | 107 |

Chapter 8: First-Order PDEs and the Method of Characteristics | 110 |

Section 8.1: The simplest PDE and the method of characteristics | 110 |

Two-dimensional graphics in MATLAB | 110 |

Section 8.2: First-order quasi-linear PDEs | 113 |

Chapter 1: Problems in multiple spatial dimensions | 14 |

Section 1.2: Fourier series on a rectangular domain | 14 |

Section 8.3: Fourier series on a disk | 117 |

Graphics on the disk | 118 |

Chapter 12: More about Fourier series | 121 |

Section 12.1: The complex Fourier series | 121 |

Section 9.2: Fourier series and the FFT | 123 |

Chapter 13: More about finite element methods | 125 |

Section 13.1 Implementation of finite element methods | 125 |

Creating a mesh | 125 |

Computing the stiffness matrix and the load vector | 128 |

Testing the code | 132 |

Introduction

MATLAB and MATLAB notebooks | I will also give a preliminary introduction to the capabilities of |

In this introduction, I will explain the organization of this tutorial and give some basic information about MATLAB.

About this tutorial

techniques presented in my textbook | This really is a tutorial (not a reference), meant to be read and used |

in parallel with the textbook | For this reason, I have structured the tutorial to have the same chapter and |

The purpose of this document is to explain the features of MATLAB that are useful for applying the sections titles as the book. However, the purpose of the sections of this document is not to re-explain the material in the text; rather, it is to present the capabilities of MATLAB as they are needed by someone studying the text.

Therefore, for example, in Section 2.1, "Heat flow in a bar; Fourier's Law", I do not explain any physics or modeling. (The physics and modeling are found in the text.) Instead, I explain the MATLAB command for integration, because Section 2.1 is the first place in the text where the student is asked to integrate a function. Because of this style of organization, some parts of the text have no counterpart in this tutorial. For example, there is no Chapter 7, because, by the time you have worked through the first six chapters of the tutorial, you have learned all of the capabilities of MATLAB that you need to address the material in Chapter 7 of the text. For the same reason, you will see that some individual sections are missing; Chapter 5, for example, begins with Section 5.2.

I should point out that my purpose is writing this tutorial is not to show you how to solve the problems in the text; rather, it is to give you the tools to solve them. Therefore, I do not give you a worked-out example of every problem type---if I did, your "studying" could degenerate to simply looking for an example, copying it, and making a few changes. At crucial points, I do provide some complete examples, since I see no other way to illustrate the power of MATLAB than in context. However, there is still plenty for you to figure out for yourself!

About MATLAB

MATLAB, which is short for Matrix Laboratory, incorporates numerical computation, symbolic computation, graphics, and programming. As the name suggests, it is particularly oriented towards matrix computations, and it provides both state-of-the-art algorithms and a simple, easy to learn interface for manipulating matrices. In this tutorial, I will touch on all of the capabilities mentioned above: numerical and symbolic computation, graphics, and programming.

This document you are reading is called an M-Book. It integrates text and MATLAB commands (with their output, including graphics). If you are running MATLAB under Microsoft Windows, then an MBook becomes an interactive document: by running the M-Book under MATLAB, you can enter new MATLAB commands and see their output inside the M-Book itself. The MATLAB command that allows you to do this is called notebook. To run this tutorial under MATLAB, just type "notebook tutorial.docx" at the MATLAB prompt. The file tutorial.docx must be in the working directory or in some directory in the MATLAB path (both of these concepts are explained below.)

Since the M-Book facility is available only under Microsoft Windows, I will not emphasize it in this tutorial. However, Windows users should take advantage of it. The most important thing to understand about a M-Book is that it is interactive---at any time you can execute a MATLAB command and see what it does. This makes a MATLAB M-Book a powerful learning environment: when you read an explanation of a MATLAB feature, you can immediately try it out.

Getting help with MATLAB commands

Documentation about MATLAB and MATLAB commands is available from within the program itself. If you know the name of the command and need more information about how it works, you can just type "help <command name>" at the MATLAB prompt. In the same way, you can get information about a group of commands with common uses by typing "help <topic name>". I will show examples of using the command-line help feature below.

The MATLAB desktop contains a help browser covering both reference and tutorial material. To access the browser, click on the Help menu and choose MATLAB Help. You can then choose "Getting Started" from the table of contents for a tutorial introduction to MATLAB, or use the index to find specific information.

Getting started with MATLAB

As mentioned above, MATLAB has many capabilities, such as the fact that one can write programs made up of MATLAB commands. The simplest way to use MATLAB, though, is as an interactive computing environment (essentially, a very fancy graphing calculator). You enter a command and MATLAB executes it and returns the result. Here is an example:

4 |

clear 2+2 ans =

You can assign values to variables for later use:

2 |

x=2 x =

The variable x can now be used in future calculations:

4 |

x2 ans =

At any time, you can list the variables that are defined with the who command: who

Your variables are: ans x

At the current time, there are 2 variables defined. One is x, which I explicitly defined above. The other is ans (short for "answer"), which automatically holds the most recent result that was not assigned to a variable (you may have noticed how ans appeared after the first command above). You can always check the value of a variable simply by typing it:

2 |

x x =

4 |

ans ans =

If you enter a variable that has not been defined, MATLAB prints an error message:

??? Undefined function or variable 'y' |

To clear a variable from the workspace, use the clear command: who

ans x |

ans |

Your variables are: clear x who Your variables are:

To clear of the variables from the workspace, just use clear by itself:

clear who

MATLAB knows the elementary mathematical functions: trigonometric functions, exponentials, logarithms, square root, and so forth. Here are some examples:

1.4142 |

sqrt(2) ans =

0.8660 |

sin(pi/3) ans =

2.7183 |

exp(1) ans = log(ans)

1 |

ans =

A couple of remarks about the above examples: MATLAB knows the number , which is called pi.

Computations in MATLAB are done in floating point arithmetic by default. For example,

MATLAB computes the sine of /3 to be (approximately) 0.8660 instead of exactly 3/2. A complete list of the elementary functions can be obtained by entering "help elfun":

help elfun Elementary math functions.

sin - Sine. | |

sind - Sine of argument in degrees. | |

sinh - Hyperbolic sine. | |

asin - Inverse sine. | |

asind - Inverse sine, result in degrees. | |

asinh - Inverse hyperbolic sine. | |

cos - Cosine. | |

cosd - Cosine of argument in degrees. | |

cosh - Hyperbolic cosine. | |

acos - Inverse cosine. | |

acosd - Inverse cosine, result in degrees. | |

acosh - Inverse hyperbolic cosine. | |

tan - Tangent. | |

tand - Tangent of argument in degrees. | |

tanh - Hyperbolic tangent. | |

atan - Inverse tangent. | |

atand - Inverse tangent, result in degrees. | |

atan2 - Four quadrant inverse tangent. | |

atanh - Inverse hyperbolic tangent. | |

sec - Secant. | |

secd - Secant of argument in degrees. | |

sech - Hyperbolic secant. | |

asec - Inverse secant. | |

asecd - Inverse secant, result in degrees. | |

asech - Inverse hyperbolic secant. | |

csc - Cosecant. | |

cscd - Cosecant of argument in degrees. | |

csch - Hyperbolic cosecant. | |

acsc - Inverse cosecant. | |

acscd - Inverse cosecant, result in degrees. | |

acsch - Inverse hyperbolic cosecant. | |

cot - Cotangent. | |

cotd - Cotangent of argument in degrees. | |

coth - Hyperbolic cotangent. | |

acot - Inverse cotangent. | |

acotd - Inverse cotangent, result in degrees. | |

acoth - Inverse hyperbolic cotangent. | |

hypot - Square root of sum of squares. |

Trigonometric.

exp - Exponential. | |

expm1 - Compute exp(x)-1 accurately. |

Exponential. log - Natural logarithm.

log1p - Compute log(1+x) accurately. | |

log10 - Common (base 10) logarithm. | |

log2 - Base 2 logarithm and dissect floating point number. | |

pow2 - Base 2 power and scale floating point number. | |

realpow - Power that will error out on complex result. | |

reallog - Natural logarithm of real number. | |

realsqrt - Square root of number greater than or equal to zero. | |

sqrt - Square root. | |

nthroot - Real n-th root of real numbers. | |

nextpow2 - Next higher power of 2. |

abs - Absolute value. | |

angle - Phase angle. | |

complex - Construct complex data from real and imaginary parts. | |

conj - Complex conjugate. | |

imag - Complex imaginary part. | |

real - Complex real part. | |

unwrap - Unwrap phase angle. | |

isreal - True for real array. | |

cplxpair - Sort numbers into complex conjugate pairs. |

Complex.

fix - Round towards zero. | |

floor - Round towards minus infinity. | |

ceil - Round towards plus infinity. | |

round - Round towards nearest integer. | |

mod - Modulus (signed remainder after division). | |

rem - Remainder after division. | |

sign - Signum. |

Rounding and remainder.

For more information about any of these elementary functions, type "help <function_name>". For a list of help topics like "elfun", just type "help". There are other commands that form part of the help system; to see them, type "help help".

MATLAB does floating point arithmetic using the IEEE standard, which means that numbers have about 16 decimal digits of precision (the actual representation is in binary, so the precision is not exactly 16 digits). However, MATLAB only displays 5 digits by default. To change the display, use the format command. For example, "format long" changes the display to 15 digits:

format long pi ans = 3.141592653589793

Other options for the format command are "format short e" (scientific notation with 5 digits) and "format long e" (scientific notation with 15 digits).

In addition to pi, other predefined variables in MATLAB include i and j, both of which represent the imaginary unit: i=j=sqrt(-1).

clear i2 ans = -1

-1 |

j2 ans =

Although it is usual, in mathematical notation, to use i and j as arbitrary indices, this can sometimes lead to errors in MATLAB because these symbols are predefined. For this reason, I will use i and j as my standard indices when needed.

Vectors and matrices in MATLAB

The default type for any variable or quantity in MATLAB is a matrix---a two-dimensional array. Scalars and vectors are regarded as special cases of matrices. A scalar is a 1 by 1matrix, while a vector is an n by 1 or 1 by n matrix. A matrix is entered by rows, with entries in a row separated by spaces or commas, and the rows separated by semicolons. The entire matrix is enclosed in square brackets. For example, I can enter a 3 by 2 matrix as follows:

1 2 | |

3 4 | |

5 6 |

A=[1 2;3 4;5 6] A =

Here is how I would enter a 2 by 1 (column) vector:

x=[1;-1]

1 | |

-1 |

x =

A scalar, as we have seen above, requires no brackets:

4 |

a=4 a =

A variation of the who command, called whos, gives more information about the defined variables:

Name | Size Bytes Class Attributes |

A | 3x2 48 double |

a | 1x1 8 double |

ans | 1x1 8 double |

x | 2x1 16 double |

whos

The column labeled "size" gives the size of each array; you should notice that, as I mentioned above, a scalar is regarded as a 1 by 1 matrix (see the entry for a, for example).

MATLAB can perform the usual matrix arithmetic. Here is a matrix-vector product:

-1 |

A*x ans = -1

-1 |

Here is a matrix-matrix product:

B=[-1 3 4 6;2 0 1 -2]

-1 3 4 6 | |

2 0 1 -2 |

3 3 6 2 | |

5 9 16 10 | |

7 15 26 18 |

A*B ans =

MATLAB signals an error if you attempt an operation that is undefined:

Inner matrix dimensions must agree |

B*A ??? Error using ==> mtimes

Matrix dimensions must agree |

A+B ??? Error using ==> plus More about M Books

If you are reading this document using the MATLAB notebook facility, then you may wish to execute the commands as you read them. Otherwise, the variables shown in the examples are not actually created in the MATLAB workspace. To execute a command, click on it (or select it) and press control-enter (that is, press the enter key while holding down the control key). While reading the tutorial, you should execute each of my commands as you come to it. Otherwise, the state of MATLAB is not what it appears to be, and if you try to experiment by entering your own commands, you might get unexpected results if your calculations depend on the ones you see in this document.

Notice that the command lines in this document appear in green, and are enclosed in gray square brackets. Output appears in blue text, also enclosed in gray square brackets. These comments may not apply if you are reading a version of this document that has been printed or converted to another format (such as or PDF).

If you are reading this using MATLAB's notebook command, then, as I mentioned above, you can try your own MATLAB commands at any time. Just move the cursor to a new line, type the command, and then type control-enter. You should take advantage of this facility, as it will make learning MATLAB much easier.

Simple graphics in MATLAB

Two-dimensional graphics are particularly easy to understand: If you define vectors x and y of equal length (each with n components, say), then MATLAB's plot command will graph the points (x1,y1), (x2,y2), …, (xn,yn) in the plane and connect them with line segments. Here is an example:

format short x=[0,0.25,0.5,0.75,1] x =

0 0.2500 0.5000 0.7500 1.0 |

1 0 1 0 1 |

y=[1,0,1,0,1] y = plot(x,y)

Two features of MATLAB make it easy to generate graphs. First of all, the command linspace creates a vector with linearly spaced components---essentially, a regular grid on an interval. (Mathematically, linspace creates a finite arithmetic sequence.) To be specific, linspace(a,b,n) creates the (row) vector whose components are a,a+h,a+ 2h,…,a+(n-1)h, where h=1/(n-1).

0 0.2000 0.4000 0.6000 0.8000 1.0 |

x=linspace(0,1,6) x =

The second feature that makes it easy to create graphs is the fact that all standard functions in MATLAB, such as sine, cosine, exp, and so forth, are vectorized. A vectorized function f, when applied to a vector x, returns a vector y (of the same size as x) with ith component equal to f(xi). Here is an example:

y=sin(pi*x)

0 0.5878 0.9511 0.9511 0.5878 0.0 |

y =

I can now plot the function: plot(x,y)

Of course, this is not a very good plot of sin( x), since the grid is too coarse. Below I create a finer grid and thereby obtain a better graph of the function. Often when I create a new vector or matrix, I do not want MATLAB to display it on the screen. (The following example is a case in point: I do not need to see the 41 components of vector x or vector y.) When a MATLAB command is followed by a semicolon, MATLAB will not display the output.

x=linspace(0,1,41); y=sin(pi*x); plot(x,y)

The basic arithmetic operations have vectorized versions in MATLAB. For example, I can multiply two vectors component-by-component using the ".*" operator. That is, z=x.*y sets zi equal to xiyi. Here is an example:

1 | |

2 | |

3 |

x=[1;2;3] x =

2 | |

4 | |

6 |

y=[2;4;6] y =

2 | |

8 | |

18 |

z=x.*y z =

The "./" operator works in an analogous fashion. There are no ".+" or ".-" operators, since addition and subtraction of vectors are defined componentwise already. However, there is a "." operator that applies an exponent to each component of a vector.

1 | |

2 | |

3 |

x =

1 | |

4 | |

9 |

x.2 ans =

Finally, scalar addition is automatically vectorized in the sense that a+x, where a is a scalar and x is a vector, adds a to every component of x. The vectorized operators make it easy to graph a function such as f(x)=x/(1+x2). Here is how it is done:

x=linspace(-5,5,41); y=x./(1+x.2); plot(x,y)

If I prefer, I can just graph the points themselves, or connect them with dashed line segments. Here are examples:

plot(x,y,'.') plot(x,y,'o') plot(x,y,'--')

The string following the vectors x and y in the plot command ('.', 'o', and '—' in the above examples) specifies the line type. it is also possible to specify the color of the graph. For more details, see "help plot".

It is not much harder to plot two curves on the same graph. For example, I plot y=x2 and y=x3 together:

x=linspace(-1,1,101); plot(x,x.2,x,x.3)

I can also give the lines different linetypes: plot(x,x.2,'-',x,x.3,'--')

Symbolic computation in MATLAB

In addition to numerical computation, MATLAB can also perform symbolic computations. However, since, by default, variables are floating point, you must explicitly indicate that a variable is intended to be symbolic. One way to do this is using the syms command, which tells MATLAB that one or more variables are symbolic. For example, the following command defines a and b to be symbolic variables:

syms a b

I can now form symbolic expressions involving these variables:

2*a*b ans = 2*a*b

Notice how the result is symbolic, not numeric as it would be if the variables were floating point variables. Also, the above calculation does not result in an error, even though a and b do not have values.

Another way to create a symbolic variable is to assign a symbolic value to a new symbol. Since numbers are, by default, floating point, it is necessary to use the sym function to tell MATLAB that a number is symbolic:

c=sym(2) c = 2

Name | Size Bytes Class Attributes |

A | 3x2 48 double |

B | 2x4 64 double |

a | 1x1 60 sym |

ans | 1x1 60 sym |

b | 1x1 60 sym |

c | 1x1 60 sym |

x | 1x101 808 double |

y | 1x41 328 double |

z | 3x1 24 double |

whos

I can do symbolic computations:

a=sqrt(c) a = 2(1/2)

You should notice the difference between the above result and the following:

1.4142 |

a=sqrt(2) a =

Name | Size Bytes Class Attributes |

A | 3x2 48 double |

B | 2x4 64 double |

a | 1x1 8 double |

ans | 1x1 60 sym |

b | 1x1 60 sym |

c | 1x1 60 sym |

x | 1x101 808 double |

whos y 1x41 328 double

z | 3x1 24 double |

Even though a was declared to be a symbolic variable, once I assign a floating point value to it, it becomes numeric. This example also emphasizes that sym must be used with literal constants if they are to interpreted as symbolic and not numeric:

a=sqrt(sym(2)) a = 2(1/2)

As a further elementary example, consider the following two commands:

sin(sym(pi)) ans = 0 sin(pi) ans = 1.2246e-016

In the first case, since π is symbolic, MATLAB notices that the result is exactly zero; in the second, both π and the result are represented in floating point, so the result is not exactly zero (the error is just roundoff error).

Using symbolic variables, I can perform algebraic manipulations.

syms x p=(x-1)*(x-2)*(x-3) p = (x - 1)*(x - 2)*(x - 3) expand(p) ans = x3 - 6*x2 + 1*x - 6 factor(ans) ans = (x - 3)*(x - 1)*(x - 2)

(Parte **1** de 4)