Partial Differential Equations in MATLAB

Partial Differential Equations in MATLAB

(Parte 2 de 4)

Integers can also be factored, though the form of the result is depends on whether the input is numeric or symbolic:

2 2 2 2 3 3

factor(144) ans = factor(sym(144)) ans = 24*32

For a numeric (integer) input, the result is a list of the factors, repeated according to multiplicityFor a

symbolic input, the output is a symbolic expression.

An important command for working with symbolic expressions is simplify, which tries to reduce an expression to a simpler one equal to the original. Here is an example:

clear syms x a b c p=(x-1)*(a*x2+b*x+c)+x2+3*x+a*x2+b*x p = 3*x + b*x + a*x2 + (x - 1)*(a*x2 + b*x + c) + x2 simplify(p) ans = 3*x - c + c*x + a*x3 + b*x2 + x2

Since the concept of simplification is not precisely defined (which is simpler, a polynomial in factored form or in expanded form a+bx+cx2+…?), MATLAB has a number of specialized simplification commands. I have already used two of them, factor and expand. Another is collect, which "gathers like terms":

p p = 3*x + b*x + a*x2 + (x - 1)*(a*x2 + b*x + c) + x2 collect(p,x) ans = a*x3 + (b + 1)*x2 + (c + 3)*x - c

By the way, the display of symbolic output can be made more mathematical using the pretty command: pretty(ans)

3 2

a x + (b + 1) x + (c + 3) x - c

Note: MATLAB's symbolic computation is based on Maple ™, a computer algebra system originally developed at the University of Waterloo, Canada, and marketed by Waterloo Maple, Inc. If you have the Extended Symbolic Math Toolbox with your installation of MATLAB, then you have access to all nongraphical Maple functions; see "help maple" for more details. However, these capabilities are not included in the standard Student Version of MATLAB, and so I will not emphasize them in this tutorial.

Manipulating functions in MATLAB

Symbolic expressions can be treated as functions. Here is an example:

syms x p=x/(1+x2) p = x/(x2 + 1)

Using the subs command, I can evaluate the function p for a given value of x. The following command substitutes 3 for every occurrence of x in p:


subs(p,x,3) ans =

The calculation is automatically vectorized:

y=linspace(-5,5,101); z=subs(p,x,y); plot(y,z)

One of the most powerful features of symbolic computation in MATLAB is that certain calculus operations, notably integration and differentiation, can be performed symbolically. These capabilities will be explained later in the tutorial.

Above I showed how to evaluate a function defined by a symbolic expression. It is also possible to explicitly define a function at the command line. Here is an example:


f=@(x)x2 f =

The function f can be evaluated in the expected fashion, and the input can be either floating point or symbolic:


f(1) ans = f(sqrt(pi))


ans = f(sqrt(sym(pi))) ans = pi

The @ operator can also be used to create a function of several variables:



g = g(1,2) ans =


20 g(pi,14) ans = 138.1745

A function defined by the @ operator is not automatically vectorized:

x=linspace(0,1,1)'; y=linspace(0,1,1)'; g(x,y) ??? Error using ==> mpower Inputs must be a scalar and a square matrix. To compute elementwise POWER, use POWER (.^) instead.

Error in ==> @(x,y)x2*y

The function can be vectorized when it is defined:


g=@(x,y)x.2.*y g = g(x,y)


ans =

There is a third way to define a function in MATLAB, namely, to write a program that evaluates the function. I will defer an explanation of this technique until Chapter 3, where I first discuss programming in MATLAB.

Saving your MATLAB session

When using MATLAB, you will frequently wish to end your session and return to it later. Using the save command, you can save all of the variables in the MATLAB workspace to a file. The variables are stored efficiently in a binary format to a file with a ".mat" extension. The next time you start MATLAB, you can load the data using the load command. See "help save" and "help load" for details.

About the rest of this tutorial

The remainder of this tutorial is organized in parallel with my textbook. Each section in the tutorial introduces any new MATLAB commands that would be useful in addressing the material in the corresponding section of the textbook. As I mentioned above, some sections of the textbook have no counterpart in the tutorial, since all of the necessary MATLAB commands have already been explained. For this reason, the tutorial is intended to be read from beginning to end, in conjunction with the textbook.

Chapter 1: Classification of differential equations

As I mentioned above, MATLAB can perform symbolic calculus on expressions. Consider the following example:

syms x f=sin(x2) f = sin(x2)

I can differentiate this expression using the diff command:

diff(f,x) ans = 2*x*cos(x2)

The same techniques work with expressions involving two or more variables:

syms x y q=x2*y3*exp(x) q = x2*y3*exp(x) pretty(q)

2 3 x y exp(x) diff(q,y) ans = 3*x2*y2*exp(x)

Thus MATLAB can compute partial derivatives just as easily as ordinary derivatives.

One use of these capabilities to test whether a certain function is a solution of a given differential equation. For example, suppose I want to check whether the function u(t)=eat is a solution of the ODE

I define syms a t u=exp(a*t) u = exp(a*t)

I can then compute the left side of the differential equation, and see if it agrees with the right side (zero):

diff(u,t)-a*u ans = 0

Thus the given function u is a solution. Is the function v(t)=at another solution? I can check it as follows:

v=a*t v = a*t diff(v,t)-a*v ans = a - a2*t

Since the result is not zero, the function v is not a solution.

It is no more difficult to check whether a function of several variables is the solution of a PDE. For example, is w(x,y)=sin(πx)+sin(πy) a solution of the differential equation

As before, I can answer this question by defining the function and substituting it into the differential equation:

syms x y w=sin(pi*x)+sin(pi*y) w = sin(pi*x) + sin(pi*y) diff(w,x,2)+diff(w,y,2) ans =

- pi2*sin(pi*x) - pi2*sin(pi*y) simplify(ans) ans = -pi2*(sin(pi*x) + sin(pi*y))

Since the result is not zero, the function w is not a solution of the PDE.

The above example shows how to compute higher derivatives of an expression. For example, here is the fifth derivative of w with respect to x:

diff(w,x,5) ans = pi5*cos(pi*x)

To compute a mixed partial derivative, we have to iterate the diff command. Here is the mixed partial derivative of w(x,y)=x2+xy2 with respect to x and then y:

syms x y w=x2*exp(y)+x*y2 w = x2*exp(y) + x*y2 diff(diff(w,x),y) ans = 2*y + 2*x*exp(y)

Instead of using expressions in the above calculations, I can use functions. Consider the following:

clear syms a x f=@(x)exp(a*x) f = @(x)exp(a*x) diff(f(x),x)-a*f(x) ans = 0

When defining a function that depends on a parameter (like the function f above, which depends on a), the value of the parametered is "captured" at the time when the function is defined.


clear syms a f=@(x)exp(a*x) f = f(1) ans = exp(a) a=1.5


a = f(1) ans = exp(a) a


The same is true if the parameter has a numeric value:


clear a=2 a = f=@(x)exp(a*x)


f = f(1)


ans = a=3


a = f(1) ans = 7.3891

Chapter 2: Models in one dimension

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

MATLAB can compute both indefinite and definite integrals. The command for computing an indefinite integral (antiderivative) is exactly analogous to that for computing a derivative:

syms x f=x2 f = x2 int(f,x) ans = x3/3

As this example shows, MATLAB does not add the "constant of integration." It simply returns one antiderivative (when possible). If the integrand is too complicated, MATLAB just returns the integral unevaluated, and prints a warning message.


Warning: Explicit integral could not be found. ans = int(exp(cos(x)), x)

To compute a definite integral, we must specify the limits of integration:

int(x2,x,0,1) ans = 1/3

MATLAB has commands for estimating the value of a definite integral that cannot be computed analytically. Consider the following example:


Warning: Explicit integral could not be found. ans = int(exp(cos(x)), x = 0..1)

Since MATLAB could not find an explicit antiderivative, I can use the quad function to estimate the definite integral. The quad command takes, as input, a function rather than an expression (as does int). Therefore, I must first create a function:


f=@(x)exp(cos(x)) f =

Now I can invoke quad:

quad(f,0,1) ans = 2.3416

"quad" is short for quadrature, another term for numerical integration. For more information about the quad command, see "help quad".

As a further example of symbolic calculus, I will use the commands for integration and differentiation to test Theorem 2.1 from the text. The theorem states that (under certain conditions) a partial derivative can be moved under an integral sign:

dcd c dyyxx FdyyxFdx

Here is a specific instance of the theorem:

syms x y c d f=x*y3+x2*y f = x2*y + x*y3 r1=diff(int(f,y,c,d),x) r1 = - (x*(c2 - d2))/2 - ((c2 - d2)*(c2 + d2 + 2*x))/4 r2=int(diff(f,x),y,c,d) r2 = -((c2 - d2)*(c2 + d2 + 4*x))/4 r1-r2 ans = ((c^2 - d^2)*(c^2 + d^2 + 4*x))/4 - ((c^2 - d^2)*(c^2 + d^2 + 2*x))/4 -

(x*(c2 - d2))/2 simplify(ans) ans = 0

Solving simple boundary value problems by integration

Consider the following BVP:

u xxdt

The solution can be found by two integrations (cf. Example 2.2 in the text). Remember, MATLAB does not add a constant of integration, so I do it explicitly:

clear syms x C1 C2 int(-(1+x),x)+C1 ans = C1 - (x + 1)2/2 int(ans,x)+C2 ans = C2 + C1*x - (x + 1)3/6 u=ans u = C2 + C1*x - (x + 1)3/6

The above function u, with the proper choice of C1 and C2, is the desired solution. To find the constants, I solve the (algebraic) equations implied by the boundary conditions. The MATLAB command solve can be used for this purpose.

The MATLAB solve command

Before completing the previous example, I will explain the solve command on some simpler examples. Suppose I wish to solve the linear equation ax+b=0 for x. I can regard this as a root-finding problem: I want the root of the function f(x)=ax+b. The solve command finds the root of a function with respect to a given independent variable:

syms f x a b f=a*x+b f = b + a*x solve(f,x) ans = -b/a

If the equation has more than one solution, solve returns the possibilities in an array:

syms f x f=x2-3*x+2; solve(f,x) ans = 1 2

As these examples show, solve is used to find solutions of equations of the form f(x)=0; only the expression f(x) is input to solve.

solve can also be used to solve a system of equations in several variables. In this case, the equations are listed first, followed by the unknowns. For example, suppose I wish to solve the equations x+y=1, 2x-y=1.

Here is the command:

syms x y s=solve(x+y-1,2*x-y-1,x,y)

x: [1x1 sym]
y: [1x1 sym]

s =

What kind of variable is the output s? If we list the variables in the workspace,

NameSize Bytes Class Attributes
C11x1 60 sym
C21x1 60 sym

whos a 1x1 60 sym

ans2x1 60 sym
b1x1 60 sym
f1x1 60 sym
s1x1 368 struct
u1x1 60 sym
x1x1 60 sym
y1x1 60 sym

we see that s is a 1 by 1 struct array, that is, an array containing a single struct. A struct is a data type with named fields that can be accessed using the syntax variable.field. The variable s has two fields:

x: [1x1 sym]
y: [1x1 sym]

s s =

The two fields hold the values of the unknowns x and y:

s.x ans = 2/3 s.y ans = 1/3

If the system has more than one solution, then the output of solve will be a struct, each of whose fields is an array containing the values of the unknowns. Here is an example: s=solve(x2+y2-1,y-x2,x,y)

x: [4x1 sym]
y: [4x1 sym]

s =

The first solution is pretty(s.x(1))

/ 1/2\1/2
| 5|
\ 2/

| ---- - 1/2 | pretty(s.y(1))

1/2 5 ---- - 1/2 2

The second solution is pretty(s.x(2))

/1/2 \1/2
| 5|

| - ---- - 1/2 | \ 2 /

You might notice that the second solution is complex (at least, the value of x is complex)This might be

pretty(s.y(2)) - ---- - 1/2 easier to see from the numerical value of the expressions, which we can see with the double command (which converts a symbolic expression to a floating point value, if possible)

0 + 1.2720i

double(s.x(2)) ans = double(s.y(2)) ans = -1.6180

Here are the remaining solutions (given in floating point) double(s.x(3)) ans = -0.7862 double(s.y(3))


ans =

0 - 1.2720i

double(s.x(4)) ans = double(s.y(4)) ans = -1.6180

If the equation to be solved cannot be solved exactly, then solve automatically tries to solve it numerically, using extended precision arithmetic (approximately 32 decimal digits):


syms x solve(x5+sym(4)*x4-x3+x2-x+1,x) ans = 4.3019656878883790704888463433324 0.40236742000277868343510597733001*sqrt(-1) + 0.49445817238673908475778249075192 - 0.67380934595293810995276180453239*sqrt(-1) - 0.34347532844254954951335931908572 0.67380934595293810995276180453239*sqrt(-1) - 0.34347532844254954951335931908572 0.49445817238673908475778249075192 - 0.40236742000277868343510597733001*sqrt(-1)

Finally, there is another way to call solve. The equation and unknowns can be given as strings, meaning that there is no need to define symbolic variables before calling solve:

solve('x2-2*x-1=0','x') ans =

2^(1/2) + 1 1 - 2(1/2)

Notice that, when specifying the equation in symbolically form, it must be expressed with the right-hand

side equal to 0, and only the left-hand side is passed to solve (that is, to solve for , we call solve as solve(f(x),x). However, when calling solve using strings instead of symbolic expressions, we can

use either form, (as above) or :

solve('x2-2*x-1','x') ans = 2^(1/2) + 1 1 - 2(1/2)

Back to the example

We can now solve for the constants of integrations in the solution to the BVP

u xxdt

Recall that the solution is of the form u u = C2 + C1*x - (x + 1)3/6

andNotice how I use the subs command to form and

We must use the boundary conditions to compute C1 and C2. The equations are

C1: [1x1 sym]
C2: [1x1 sym]

s=solve(subs(u,x,0)-2,subs(u,x,1),C1,C2) s =

Here are the values of C1 and C2:

s.C1 ans = -5/6 s.C2 ans = 13/6

Here is the solution:

u=subs(u,{C1,C2},{s.C1,s.C2}) u = 13/6 - (x + 1)3/6 - (5*x)/6

Notice how, when substituting for two or more variables, the variables and the values are enclosed in curly braces.

Let us check our solution:

-diff(u,x,2) ans = x + 1 subs(u,x,sym(0)) ans = 2 subs(u,x,sym(1)) ans = 0

The differential equation and the boundary conditions are satisfied. Another example

Now consider the BVP with a nonconstant coefficient:

xdx duxdxd

Integrating once yields

(It is easier to perform this calculation in my head than to ask MATLAB to integrate 0.) I now perform the second integration:

clear syms C1 C2 x u=int(C1/(1+x/2),x)+C2 u = C2 + 2*C1*log(x + 2)

Now I use solve to find C1 and C2:

C1: [1x1 sym]
C2: [1x1 sym]

s=solve(subs(u,x,0)-20,subs(u,x,1)-25,C1,C2) s =

Here is the solution:

u=subs(u,{C1,C2},{s.C1,s.C2}) u = (45035996273704960*log(x + 2))/3652105019575333 + 41825526550679865/3652105019575333

Notice the unexpected answer; the problem is that MATLAB did not interpret the constants in the equations (0, 1, 20, 25) as symbolic quantities, but rather as numerical (floating point) values. As a result,

it used extended precision arithmetic while finding a solutionThe desired solution can be found by

31 specifying that the various numbers be treated as symbolic:

clear u syms C1 C2 x u=int(C1/(1+x/2),x)+C2 u = C2 + 2*C1*log(x + 2)

C1: [1x1 sym]
C2: [1x1 sym]

s=solve(subs(u,x,sym(0))-sym(20),subs(u,x,sym(1))-sym(25),C1,C2) s = u=subs(u,{C1,C2},{s.C1,s.C2}) u = (5*(5*log(2) - 4*log(3)))/(log(2) - log(3)) - (5*log(x + 2))/(log(2) - log(3))

Now I will check the answer:

simplify(-diff((1+x/2)*diff(u,x),x)) ans = 0


subs(u,x,0) ans = subs(u,x,1) ans = 25.0

Chapter 3: Essential linear algebra

Section 3.1 Linear systems as linear operator equations

I have already showed you how to enter matrices and vectors in MATLAB. I will now introduce a few more elementary operations on matrices and vectors, and explain how to extract components from a vector and entries, rows, or columns from a matrix. At the end of this section, I will describe how MATLAB can perform symbolic linear algebra; until then, the examples will use floating point arithmetic.

(Parte 2 de 4)