Partial Differential Equations in MATLAB

Partial Differential Equations in MATLAB

(Parte 4 de 4)

Here is a sample script. Its purpose is to graph the function

on the interval [0,1]. (Recall that MATLAB cannot compute the integral explicitly, so this is a nontrivial task.) (Caveat: I did not try to make the following script efficient; indeed, it is decidedly inefficient!) type scriptex

% Define the integrand g=@(s)exp(cos(s)); % Create a grid x=linspace(0,1,21); y=zeros(1,21);

% Evaluate the integral int(exp(cos(s)),s,0,x) for % each value of x on the grid:

for i=1:length(x) y(i)=quad(g,0,x(i)); end

% Now plot the result plot(x,y) Here is the result of running scriptex.m:

clear scriptex

As I mentioned above, scripts do not have local variables. Any variables defined in scriptex exist in the MATLAB workspace:

NameSize Bytes Class Attributes
g1x1 16 function_handle
i1x1 8 double
x1x21 168 double
y1x21 168 double

whos

Section 3.4: Orthogonal bases and projection Consider the following three vectors:

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

These vectors are orthonormal, as can easily be checked. Therefore, I can easily express any vector as a linear combination of the three vectors, which form an orthonormal basis. To test this, I will use the randn command to generate a random vector (for more information, see "help randn"):

x=randn(3,1) x = -0.2050

1.4897

-0.1241 y=(x'*v1)*v1+(x'*v2)*v2+(x'*v3)*v3 y = -0.2050

1.4897

-0.1241 y-x

0.4718
0

ans = 1.0e-015 * -0.0139

Notice that the difference between y and x (which should be equal) is due to roundoff error, and is very small.

Working with the L2 inner product

Since MATLAB can compute integrals symbolically, we can use it to compute the L2 inner product and norm. For example:

clear syms x f=x f = x g=x2 g = x2 int(f*g,x,0,1) ans = 1/4

the L2 inner productThe M-file l2ip.m does this:

If you are going to perform such calculations repeatedly, it is convenient to define a function to compute help l2ip I=l2ip(f,g,a,b,x)

This function computes the L^2 inner product of two functions
f(x) and g(x), that is, it computes the integral from a to b of
f(x)*g(x). The two functions must be defined by symbolic
expressions f and g.
The variable of integration is assumed to be x. A different
variable can passed in as the (optional) fifth input.
The inputs a and b, defining the interval [a,b] of integration,
are optional. The default values are a=0 and b=1.

Notice that I assigned the default interval to [0,1]. Here is an example:

syms x l2ip(x,x2) ans = 1/4

For convenience, I also define a function computing the L2 norm:

help l2norm I=l2norm(f,a,b,x)

This function computes the L^2 inner product of the function f(x)
The functions must be defined by the symbolic expressions f.
The variable of integration is assumed to be x. A different
variable can passed in as the (optional) fourth input.
The inputs a and b, defining the interval [a,b] of integration,
are optional. The default values are a=0 and b=1.

l2norm(x) ans = 3(1/2)/3 double(ans) ans =

0.5774

Now consider the following two functions:

clear syms x pi f=x*(1-x) f = -x*(x - 1) g=8/pi3*sin(pi*x) g = (8*sin(pi*x))/pi3

(I must tell MATLAB that pi is to be regarded as symbolic.) The following graph shows that the two functions are quite similar on the interval [0,1]:

t=linspace(0,1,51); f1=subs(f,x,t); g1=subs(g,x,t); plot(t,f1,'-',t,g1,'--')

By how much do the two functions differ? One way to answer this question is to compute the relative difference in the L2 norm:

l2norm(f-g)/l2norm(f) ans = 30(1/2)*(1/30 - 32/pi6)(1/2) double(ans) ans = 0.0380

The difference is less than 4%.

Here are two more examples from Section 3.4 that illustrate some of the capabilities of MATLAB. Example 3.37

points (xi,yi)The data given in this example can be stored in two vectors:

The purpose of this example to compute the first-degree polynomial f(x)=mx+b best fitting given data clear x=[0.1;0.3;0.4;0.75;0.9]; y=[1.7805;2.2285;2.3941;3.2226;3.5697];

Here is a plot of the data: plot(x,y,'o')

To compute the first-order polynomial f(x)=mx+b that best fits this data, I first form the Gram matrix. The ones command creates a vector of all ones:

1
1
1
1
1

e=ones(5,1) e =

1.6325 2.4500

G=[x'*x,x'*e;e'*x,e'*e] G = 2.4500 5.0

Next, I compute the right hand side of the normal equations:

7.4339

z=[x'*y;e'*y] z = 13.1954

Now I can solve for the coefficients c=[m;b]:

2.2411
1.5409

c=G\z c = I will define the solution as a function:

@(x)c(1)*x+c(2)

l=@(x)c(1)*x+c(2) l =

Here is a plot of the best fit line, together with the data:

t=linspace(0,1,1); plot(t,l(t),x,y,'o')

The fit is not bad. Example 3.38

In this example, I will compute the best quadratic approximation to the function ex on the interval [0,1]. Here are the basis functions for the subspace P2:

(Parte 4 de 4)

Comentários