Section 3.4: Orthogonal bases and projection

Working with the L² inner product

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

syms x

f=x

f =

g =

ans =

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.

l2ip(x,x^2)

ans =

1/4
For convenience, I also define a function computing the L² norm:

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.

ans =

double(ans)

ans =

0.5774

Example 3.35

syms x pi

f=x*(1-x)

f =

g =

f1=subs(f,x,t);

g1=subs(g,x,t);

plot(t,f1,'-',t,g1,'--')

ans =

double(ans)

ans =

0.0380
The difference is less than 4%.

Example 3.37

The purpose of this example to compute the first-degree polynomial f(x)=mx+b best fitting given data points (x_i,y_i). The data given in this example can be stored in two vectors:

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:
e=ones(5,1)

e =

1

1

1

1

1
G=[x'*x,x'*e;e'*x,e'*e]

G =

1.6325 2.4500

2.4500 5.0000
Next, I compute the right hand side of the normal equations:
z=[x'*y;e'*y]

z =

7.4339

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

c=G\z

c =

2.2411

1.5409
I will define the solution as a function:

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

l =

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

Here is a plot of the best fit line, together with the data:
t=linspace(0,1,11);

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 e^x on the interval [0,1]. Here are the basis functions for the subspace P₂:

clear

syms x

p1=1;

p2=x;

p3=x^2;
I now compute the Gram matrix and the right hand side of the normal equations:
G=[l2ip(p1,p1),l2ip(p1,p2),l2ip(p1,p3)

l2ip(p2,p1),l2ip(p2,p2),l2ip(p2,p3)

l2ip(p3,p1),l2ip(p3,p2),l2ip(p3,p3)]

G =

[ 1, 1/2, 1/3]

[ 1/2, 1/3, 1/4]

[ 1/3, 1/4, 1/5]

b=[l2ip(p1,exp(x));l2ip(p2,exp(x));l2ip(p3,exp(x))]

b =

exp(1) - 1

1

exp(1) - 2

Now I solve the normal equations and find the best fit quadratic:
c=G\b

c =

39*exp(1) - 105

588 - 216*exp(1)

210*exp(1) - 570
Here is the solution:
q=c(1)*p1+c(2)*p2+c(3)*p3

q =

(210*exp(1) - 570)*x^2 + (588 - 216*exp(1))*x + 39*exp(1) - 105
Here is the graph of y=e^x and the quadratic approximation:
t=linspace(0,1,41)';

plot(t,exp(t),'-',t,subs(q,x,t),'--')

Since the approximation is so accurate, it is more informative to graph the error:
plot(t,exp(t)-subs(q,x,t))


We can also judge the fit by computing the relative error in the L² norm:
l2norm(exp(x)-q)/l2norm(exp(x))

ans =

(1350*exp(1) - (497*exp(2))/2 - 3667/2)^(1/2)/(exp(2)/2 - 1/2)^(1/2)

double(ans)

ans =

0.0030
The error is less than 0.3%.

Download 2.45 Mb.

Share with your friends: