Partial Differential Equations in MATLAB

(Parte 3 de 4)

clear

Consider the matrix

 1 2 3 4 5 6 7 8 9

A=[1 2 3;4 5 6;7 8 9] A =

The transpose is indicated by a single quote following the matrix name:

A' ans =

 1 4 7 2 5 8 3 6 9

Recall that, if x and y are column vectors, then the dot product of x and y can be computed as xTy:

 1 0 -1

x=[1;0;-1] x =

 1 2 3

y=[1;2;3] y =

 -2

x'*y ans =

Alternatively, I can use the dot function:

 -2

dot(x,y) ans =

I can extract a component of a vector,

 1

x(1) ans = or an entry of a matrix:

A(2,3)

 6

ans =

In place of an index, I can use a colon, which represents the entire range. For example, A(:,1) indicates all of the rows in the first column of A. This yields a column vector:

A(:,1)

 1 4 7

ans =

Similarly, I can extract a row:

A(2,:)

 4 5 6

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

MATLAB can find a basis for the null space of a matrix. Consider the matrix

 1 2 3 4 5 6 7 8 9

B=[1 2 3;4 5 6;7 8 9] B =

Here is a basis for the null space:

 0.8165

x=null(B) x = -0.4082 -0.4082

Since MATLAB returned a single vector, this indicates that the null space is one-dimensional. Here is a check of the result:

B*x ans = 1.0e-015 * -0.2220

 0

-0.4441

Notice that, since the computation was done in floating point, the product is not exactly zero, but just very close.

If a matrix is nonsingular, its null space is trivial:

 1 2 2 1

A=[1,2;2,1] A = null(A) ans = Empty matrix: 2-by-0

On the other hand, the null space can have dimension greater than one:

A=[1 1 1;2 2 2;3 3 3;]

 1 1 1 2 2 2 3 3 3
 -0.8165 0 0.4082 0.7071 0.4082 -0.7071

null(A) ans = The matrix A has a two-dimensional null space.

MATLAB can compute the inverse of a nonsingular matrix:

 1 0 -1 3 2 1 2 -1 1

A=[1 0 -1;3 2 1;2 -1 1] A =

The command is called inv:

 0.3000 0.1000 0.2000 -0.1000 0.3000 -0.4000 -0.7000 0.1000 0.2000

Ainv=inv(A) Ainv =

 1.0 -0.0 0.0 0 1.0 0.0 0 0 1.0

A*Ainv ans =

Using the inverse, you can solve a linear system

 1 1 1

b=[1;1;1] b =

 0.6

x=Ainv*b x = -0.2000

-0.4000

A*x

 1 1 1

ans =

On the other hand, computing and using the inverse of a matrix A is not the most efficient way to solve Ax=b. It is preferable to solve the system directly using some variant of Gaussian elimination. The backslash operator indicates to MATLAB that a linear system is to be solved:

 0.6

x1=A\b x1 = -0.2000

-0.4000 x-x1

 0

ans = 0

 0
 multiplying on the left by However, MATLAB does not compute the inverse.)

(To remember how the backslash operator works, just think of A\b as "dividing on the left by ", or

Section 3.3: Basis and dimension

In this section, I will further demonstrate some of the capabilities of MATLAB by repeating some of the examples from Section 3.3 of the text.

Example 3.25

Consider the three vectors v1, v2, v3 defined as follows:

 0.5774 0.5774 0.5774

v1=[1/sqrt(3);1/sqrt(3);1/sqrt(3)] v1 = v2=[1/sqrt(2);0;-1/sqrt(2)]

 0.7071 0

v2 = -0.7071 v3=[1/sqrt(6);-2/sqrt(6);1/sqrt(6)]

 0.4082
 0.4082

v3 = -0.8165

I will verify that these vectors are orthogonal:

 0

v1'*v2 ans = v1'*v3

 0

ans = v2'*v3

 0

ans = Example 3.27

I will verify that the following three quadratic polynomials for a basis for P2. Note that while I did the previous example in floating point arithmetic, this examples requires symbolic computation.

 1

clear syms p1 p2 p3 x p1=1 p1 = p2=x-1/2 p2 = x - 1/2 p3=x2-x+1/6 p3 = x2 - x + 1/6

Now suppose that q(x) is an arbitrary quadratic polynomial:

syms q a b c q=a*x2+b*x+c q = a*x2 + b*x + c

I want to show that q can be written in terms of p1, p2, and p3: syms c1 c2 c3 q-(c1*p1+c2*p2+c3*p3) ans = c - c1 + b*x + a*x2 - c2*(x - 1/2) - c3*(x2 - x + 1/6)

I need to gather like terms in this expression, which is accomplished with the collect command:

collect(ans,x) ans = (a - c3)*x2 + (b - c2 + c3)*x + c - c1 + c2/2 - c3/6

I can now set each coefficient equal to zero and solve for the coefficients:

r=solve(a-c3,b-c2+c3,c-c1+c2/2-c3/6,c1,c2,c3)

 c1: [1x1 sym] c2: [1x1 sym] c3: [1x1 sym]

r = r.c1 ans = a/3 + b/2 + c r.c2 ans = a + b r.c3 ans = a

Check:

q-(r.c1*p1+r.c2*p2+r.c3*p3) ans = b*x - b/2 - a/3 - (a + b)*(x - 1/2) + a*x2 - a*(x2 - x + 1/6) simplify(ans) ans = 0

The fact that there is a unique solution to this system, regardless of the values of a, b, c, shows that every quadratic polynomial can be uniquely written as a linear combination of p1, p2, p3, and hence that these three polynomials form a basis for P2.

Example

Here is a final example. Consider the following three vectors in R3:

u1=[1;0;2]; u2=[0;1;1]; u3=[1;2;-1];

 x=[8;2;-4];

I will verify that {u1,u2,u3} is a basis for , and express the vector in terms of the basis. As discussed in the text, {u1,u2,u3} is a basis if and only if the matrix A whose columns are u1, u2, u3 is nonsingular. It is easy to form the matrix A:

 1 0 1 0 1 2 2 1 -1

A=[u1,u2,u3] A =

One way (that works fine for such a small matrix, but is not a good method in general) to determine if A is nonsingular is to compute its determinant:

 -5

det(A) ans =

Another way to determine whether A is nonsingular is to simply try to solve a system involving A, since MATLAB will print a warning or error message if A is singular:

 3.6
 4.4

c=A\x c = -6.8000

Here is a verification of the results: x-(c(1)*u1+c(2)*u2+c(3)*u3)

 0 0 0

ans =

Symbolic linear algebra

Recall that MATLAB performs its calculations in floating point arithmetic by default. However, if desired, we can perform them symbolically. For an illustration, I will repeat the previous example.

clear syms u1 u2 u3 u1=sym([1;0;2]); u2=sym([0;1;1]); u3=sym([1;2;-1]); A=[u1,u2,u3]

A = [ 1, 0, 1]

[ 0, 1, 2]

[ 2, 1, -1] x=sym([8;2;-4]); c=A\x c = 18/5 -34/5 2/5

The solution is the same as before. Programming in MATLAB, part I

Before I continue on to Section 3.4, I want to explain simple programming in MATLAB---specifically, how to define new MATLAB functions. I have already shown you one way to do this: using the @ operator. (Also, a symbolic expression can be used in place of a function for many purposes. For instance, it can be evalutated using the subs command.) However, in addition to the @ operator, there is another method for defining MATLAB functions.

Defining a MATLAB function in an M file.

An M-file is a plain text file whose name ends in ".m" and which contains MATLAB commands. There are two types of M-files, scripts and functions. I will explain scripts later. A function is a MATLAB subprogram: it accepts inputs and computes outputs using local variables. The first line in a function must be of the form function [output1,output2,…]=function_name(input1,input2,…)

If there is a single output, the square brackets can be omitted. Also, a function can have zero inputs and/or zero outputs.

 Here is the simplest type of example Suppose I wish to define the function f(x)=sin(x2). The following

lines, stored in the M-file f.m, accomplish this:

function y=f(x) y=sin(x.2);

(Notice how I use the "." operator to vectorize the computation. Predefined MATLAB functions are always vectorized, and user-defined functions typically should be as well.) I can now use f just as a predefined function such as sin. For example:

clear x=linspace(-3,3,101); plot(x,f(x))

A few important points:

The names of user-defined functions can be listed using the what, command (similar to who, but instead of listing variables, it lists M-files in the working directory).

The contents of an M-file can be displayed using the type command.

In order for you to use an M-file, MATLAB must be able to find it, which means that the M-file must be in a directory on the MATLAB path. The MATLAB path always includes the working directory, which can be determined using the pwd (print working directory) command. The MATLAB path can be listed using the path command. Other directories can be added to the MATLAB path using the addpath command. For more information, see "help addpath".

The current path is path fempack C:\Users\Guest\Documents\MATLAB

C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\general C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\ops C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\lang C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\elmat C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\randfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\elfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\specfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\matfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\datafun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\polyfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\funfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\sparfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\scribe C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\graph2d C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\graph3d

C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\specgraph C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\graphics C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\uitools C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\strfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\imagesci C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\iofun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\audiovideo C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\timefun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\datatypes C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\verctrl C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\codetools C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\helptools C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\winfun C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\winfun\NET C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\demos C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\timeseries C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\hds C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\guide C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\plottools C:\Program Files (x86)\MATLAB\R2011a\toolbox\local C:\Program Files (x86)\MATLAB\R2011a\toolbox\matlab\datamanager C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\simulink C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\instrument C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\asynciolib C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\comparisons

C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\controllib\general C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\controllib\graphics

C:\Program Files (x86)\MATLAB\R2011a\toolbox\optim\optim C:\Program Files (x86)\MATLAB\R2011a\toolbox\optim\optimdemos C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\optimlib C:\Program Files (x86)\MATLAB\R2011a\toolbox\pde C:\Program Files (x86)\MATLAB\R2011a\toolbox\pde\pdedemos C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\eml\eml C:\Program Files (x86)\MATLAB\R2011a\toolbox\symbolic\symbolic C:\Program Files (x86)\MATLAB\R2011a\toolbox\symbolic\symbolicdemos

C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\testmeaslib\general C:\Program Files (x86)\MATLAB\R2011a\toolbox\shared\testmeaslib\graphics

The working directory is pwd ans = C:\Users\Guest\Documents\MATLAB\msgocken

The files associated with this tutorial are (on my computer) in C:\Users\Guest\Documents\MATLAB\msgocken. (When you install the tutorial on your computer, you will have to make sure that all of the files that come with the tutorial are in an accessible directory.) Here are all of the M-files in the directory called msgocken:

what msgocken

M-files in directory C:\Users\Guest\Documents\MATLAB\msgocken

 HWProblem6Function f1 l2ip scriptex beuler f2 l2norm startup eip f2a mkpl testfun euler f6 myplot testit euler1 g myplot1 testsym euler2 g1 mysubs f h mysubs1

I can look at f.m: type f function y=f(x) y=sin(x.2);

One important feature of functions defined in M-files is that variables defined in an M-file are local; that is, they exist only while the M-file is being executed. Moreover, variables defined in the MATLAB workspace (that is, the variables listed when you give the who command at the MATLAB prompt) are not accessible inside of an M-file function. Here is are examples:

type g function y=g(x)

 0.1411

a=3; y=sin(a*x); clear g(1) ans =

 ??? Undefined function or variable 'a'

The variable a is not defined in the MATLAB workspace after g executes, since a was local to the M-file g.m.

On the other hand, consider: type h function y=g(x) y=sin(a*x); clear a=3 a = 3 h(1) ??? Undefined function or variable 'a'.

Error in ==> h at 3 y=sin(a*x);

The M-file h.m cannot use the variable a, since the MATLAB workspace is not accessible within h.m. (In short, we say that a is not "in scope".)

Here is an example with two outputs: type f1 function [y,dy]=f1(x) y=3*x.^2-x+4; dy=6*x-1;

The M-file function f1.m computes both f(x)=3x2-x+4 and its derivative. It can be used as follows:

 6
 5

[v1,v2]=f1(1) v1 = v2 =

As another example, here is a function of two variables: type g1 function z=g1(x,y) z=2*x.2+y.2+x.*y;

 8

g1(1,2) ans =

A function that is defined in an m-file can be given an alias---another name---inside of MATLAB. This is done by using the @ operator to create a "function handle". This is useful because it allows you to give the file a meaningful name, such as HWproblem6Function.m, and then use a convenient alias, such as f, when typing MATLAB commands.

type HWProblem6Function function y=HWProblem6Function(x) y=exp(cos(x)).*sin(x).2 - exp(cos(x)).*cos(x);

 @HWProblem6Function

f=@HWProblem6Function f = x=linspace(-6*pi,6*pi,401)'; y=f(x); plot(x,y)

As we will see later, function handles are useful for passing one function as input to another function. Optional inputs with default values

I now want to explain the use of optional inputs in M-files. Suppose you are going to be working with the function f2(x)=sin(ax2), and you know that, most of the time, the parameter a will have value 1. However, a will occasionally take other values. It would be nice if you only had to give the input a when its value is not the typical value of 1. The following M-file has this feature:

type f2 function y=f2(x,a) if nargin<2 a=1; end y=sin(a.*x.2);

The first executable statement, "if nargin<2", checks to see if f2 was invoked with one input or two. The keyword nargin is an automatic (local) variable whose value is the number of inputs. The M-file f2.m assigns the input a to have the value 1 if the user did not provide it. Now f2 can be used with one or two parameters:

f2(pi) ans = -0.4303 f2(pi,sqrt(2)) ans =

 0.9839

If the percent sign % is used in a MATLAB command, the remainder of the line is considered to be a comment and is ignored by MATLAB. Here is an example:

if sin(pi)==0 % Testing for roundoff error disp('No roundoff error') else disp('Roundoff error detected') end

(Notice the use of the disp command, which displays a string or the value of a variable. See "help disp" for more information. Notice also the use of the if-else block, which I discuss later in the tutorial.)

The most common use of comments is for documentation in M-files. Here is a second version of the function f2 defined above:

type f2a function y=f2a(x,a)

%y=f2a(x,a) %

% This function implements the function f2(x)=sin(a*x). The parameter

% a is optional; if it is not provided, it is taken to be 1.

if nargin<2 a=1; end y=sin(a*x);

Notice how the block of comment lines explains the purpose and usage of the function. One of the convenient features of the MATLAB help system is that the first block of comments in an M-file is displayed when "help" is requested for that function:

help f2a y=f2a(x,a)

 a is optional; if it is not provided, it is taken to be 1

This function implements the function f2(x)=sin(a*x). The parameter

I will explain more about MATLAB programming in Chapter 4. M files as scripts

 MATLAB prompt. Scripts do not have local variables, and do not accept inputs or return outputs A

An M-file that does not begin with the word function is regarded as a script, which is nothing more than a collection of MATLAB commands that are executed sequentially, just as if they had been typed at the common use for a script is to collect the commands that define and solve a certain problem (e.g. a homework problem).

(Parte 3 de 4)