Least Squares with Orthogonal basis

I'm asked to do the following with p(t) = a+ bt +ct^2 + dt^3 First, determine by hand the orthogonal basis of P^3 on the interval [−1, 1] by using of the Gram-Schmidt procedure. Then transform the derived basis to the interval [1900,1990] by the simple change of variables x = (t − 1945)/45. Use this basis to: (a)Determine the rectangular matrix A and right-hand side b of the least squares problem. (b) Form and solve the Normal Equations. Use the commands A’ to transpose A and \ to solve the linear system. Plot the table and p(x) using plot(t,y,’.’,x,p(x)), where x=[1900:0.01:1990].

I have that

%p0 = 1 
%p1 = (t-1945)/45 
%p2 = ((t-1945))/45^2 - 1/3 
%p3 = ((t-1945)/45)^3 - 3/5((t-1945)/45) 

t = [1900; 1910; 1920; 1930; 1940; 1950; 1960; 1970; 1980; 1990]; 
b = [75.995; 91.972; 105.711; 123.203; 131.669; 150.697; 179.323; 203.212; 226.505; 249.633]; 
x = (1900:0.01:1990); 

r = (t-1945)/2; 
r2 = r.^2 - 1/3; 
r3 = r.^3 -(3/5)*r; 
O = ones(10,1); 
A = [O, r, r2 r3 ]; 
K = A'*A; 
Ab = A'*b; 
v = K\Ab; 

q(x) = v(1)+v(2)*((x-1945)/(45))+v(3)*[((x-1945)/(45)).^2-(1/3)]+v(4)*[((x-1945)/(45)).^3-(3/5)*((x-1945)/(45))] 
plot(t,b,'.',x,q)

My trouble is evaluating q using the coefficients v and plotting it, it's yelling at me saying "Subscript indices must either be real positive integers or logicals." and if I try to do. For the non-orthogonal basis I was just using polyval, but I wasn't sure how to apply it to the orthogonal case.

Answer 1

[John D'Errico]

Oh, this problem again. I remember you. lol. But at least you are making a credible effort on your work, so well done.

Why are you trying to define the variable as q(x)?

See that x is a vector that does not contain integers, in fact, very non-integer numbers.

q(x) attempts to index/subscript the variable q with the vector x. So of course it must fail, because as the error message says:

"Subscript indices must either be real positive integers or logicals."

However, why use a subscript at all? You want to create the vector q. So create it! (Yes, you THINK of q as a function of x. However, that is not how you would define q as a variable.)

q = v(1)+v(2)*((x-1945)/(45))+v(3)*[((x-1945)/(45)).^2-(⅓)]+v(4)*[((x-1945)/(45)).^3-(3/5)*((x-1945)/(45))];

As an aside, I would actually create the functions p0, p1, p2, p3. This will make your code better, easier to use, because you need not type those functions over and over again.

p0 = @(t) ones(size(t));
p1 = @(t) (t-1945)/45;
p2 = @(t) ((t-1945)/45).^2 - 1/3;
p3 = @(t) ((t-1945)/45).^3 - 3/5*((t-1945)/45);

(Note that I corrected a typo in how you had written p2 in your comment. Oh, you also have a typo in the expression for r, where you divide by 2, not by 45. You should proofread your code more carefully.)

Now you can write q simply as:

q = v(1)*p0(x) + v(2)*p1(x) + v(3)*p2(x) + v(4)*p3(x);

and create A as:

A = [p0(t),p1(t),p2(t),p3(t)];

Finally, it is a TERRIBLE idea to use the variable O. See how easily it can be confused with the number 0. Now imagine how a typo will one day cause you an immeasurable amount of pain when you try to find that typing error. And since O and 0 are so close on a standard keyboard, this is a serious risk.

[Megan]

Thank you so much for the help and for catching my typos, I got everything working and I'm finally done! Time to focus on finals!

[John D'Errico]

Good luck on them. I'm always happy to help when someone makes a serious effort.