$Post:=If[MatrixQ[#],MatrixForm[#],#]&

Problems on Steady State Vectors

Problem 1
a)

M={{2/3,4/9},{1/3,5/9}}

s[0]={1,0}

Do[                   \
   s[i]=M.s[i-1];     \
   Print[s[i]//N],    \
   {i,5}              \
  ]

b)

s[0]={0,1}

Do[                   \
   s[i]=M.s[i-1];     \
   Print[s[i]//N],    \
   {i,5}              \
  ]

Problem 2

M=1/10{{2,3,2,4,2},   \
       {3,2,1,1,2},   \
       {1,1,3,1,2},   \
       {3,2,3,1,3},   \
       {1,2,1,3,1}}

s[0]=Array[1/5&,5]

Do[                   \
   s[i]=M.s[i-1];     \
   Print[s[i]//N],    \
   {i,10}             \
  ]

In 9 "years", the steady state vector stabilizes
over 6 digits of precision.

Problem 3

A=Array[Max,{10,10}]

d=Table[Sum[A[[i,j]],{j,10}],{i,10}]

M=Table[A[[i,j]]/d[[j]],{i,10},{j,10}]

NullSpace[IdentityMatrix[10]-M]

s=%/Sum[%[[1,j]],{j,10}]

Problems on Orthogonal Projection

<<LinearAlgebra`Master`

Problem 1

A=Table[1/(i+j),{i,7},{j,5}]

v=Table[i-1,{i,7}]

RowReduce[Transpose[A]]

A has linearly independent columns.

GramSchmidt[Transpose[A],Normalized->False]

Reason why we do not normalize is an important amount
of time being saved, especially in projecting and
taking the sums. Results would be the same.

u=Sum[Projection[v,%[[i]]],{i,5}]

LinearSolve[A,u]

Indeed, we can check that Ax=v has no solution itself:

LinearSolve[A,v]

Problem 2

A=Array[Min,{5,7}]

v=Table[i,{i,5}]

NullSpace[A]

sh=%[[1]]t1+%[[2]]t2

sp=LinearSolve[A,v]

sg=sh+sp

u=A.sg//Together

Check:

RowReduce[Transpose[A]]

u=Sum[Projection[v,%[[i]]],{i,5}]

We see that the projection of v is itself, since v is
actually in the column space of A.

Problem 3

A=Array[#1+#2&,{5,5}]

v=Array[#&,5]

a)

NullSpace[A]

sh=%[[1]]t1+%[[2]]t2+%[[3]]t3

sp=LinearSolve[A,v]

u=sh+sp

c=A.u//Together

Again here, c=v.

b)

NullSpace[Transpose[A]]

As an observation, in this particular case, A=Transpose[A].

sh=%[[1]]s1+%[[2]]s2+%[[3]]s3

sp=LinearSolve[Transpose[A],u/.{t1->-2/5}/.{t2->-1/10}/.{t3->1/5}]

Note that there is actually only one set of t1,t2,t3 for
which the above system has any solution at all, but
Mathematica does not come up with it by itself.

sg=sh+sp

r=Transpose[A].sg//Together

c)

A.r-c

d)

In order to apply the hint, we denote by w a solution
of Transpose[A].x=c.
Any w can be written as the sum of r (r in the row
space of A, the orthogonal projection of w) and a vector
rr, orthogonal to the row space of A. Since r and rr are
orthogonal, we have ||w||^2=||r||^2+||rr||^2 ==> 
||w|| >= ||r|| and since r is unique, it is indeed the
solution of least norm.


Problem 4

Obviously Transpose[A].A is (n*m).(m*n) = n*n.
Ax=0 ==> Transpose[A].Ax=0
Transpose[A].Ax=0 ==> 0=(Transpose[A].Ax).x=(Ax).(Ax) ==> Ax=0
Therefore Ax=0 <==> Transpose[A].Ax=0 and
rank[A]=n ==> rank[Transpose[A].A]=n ==>
Transpose[A].A invertible.

A=Array[Min,{5,4}]

b=Array[#^2&,5]

x0=Inverse[Transpose[A].A].Transpose[A].b

Problems on Curve Fitting

Problem 1

x=Table[2 Pi i/20,{i,20}]

y=Table[Sin[x[[i]]],{i,20}]

M=Table[x[[i]]^(j-1),{i,20},{j,6}]; M//TableForm

B=GramSchmidt[Transpose[M],Normalized->False];

b=Sum[Projection[y//N,B[[i]]//N],{i,6}]

The reason to work with floating point representations
is that the symbolic calculation with large matrices
requires more memory than is usually available on
personal computers, whereis the use of virtual memory
takes an unreasonable amount of run-time. Since we are 
looking for an approximating polynomial (nevertheless
optimal), this should not bother a great deal.

a=LinearSolve[Array[M[[#]]&,6]//N,Array[b[[#]]&,6]]

Keeping only the first 6 equations of the system is 
necessary since, because of the truncation errors of
floating point representations, the redundant system
becomes actually unsolvable.

p[t_]:=Sum[a[[i]]t^(i-1),{i,6}]; Clear[t]; p[t]

lp=ListPlot[Table[{x[[i]],y[[i]]},{i,20}]]

graph=Plot[p[t],{t,Pi/10,2Pi}]

Show[graph,lp]

We can actually see that the fit is actually much better
than the degree-3 fit presented in the text example.

Problem 2

x=Array[#&,10]

q[u_]:=Product[i(u-3),{i,5}]; Clear[u]; q[u]//Expand

y=Table[q[x[[i]]],{i,10}]

M=Table[x[[i]]^(j-1),{i,10},{j,4}]

GramSchmidt[Transpose[M],Normalized->False];

b=Sum[Projection[y,%[[i]]],{i,4}]//Together

a=LinearSolve[M,b]//Together

p[t_]:=Sum[a[[i]]t^(i-1),{i,4}]; Clear[t]; p[t]

lp=ListPlot[Table[{x[[i]],y[[i]]},{i,10}]]

graph=Plot[p[t],{t,1,10},PlotRange->All]

We use PlotRange->All in order to tell Mathematica that
we don't want a plot with default optimized range, but
rather a complete one, in order to be able to compare
it with the previous one (the single points).

Show[graph,lp]

Min[Abs[p[x]-y],{i,10}]

Max[Abs[p[x]-y],{i,10}]

Problems on Curve Fitting with Other Inner Products

Problem 1

Clear[x,v,w]
Do[v[i]=x^(i-1),{i,8}]
iprod[u_,v_]:=Integrate[u*v,{x,0,Pi}]
Projection[u_,v_]:=iprod[u,v]/iprod[v,v]v
Do[w[i]=v[i]-Sum[Projection[v[i],w[j]],{j,i-1}],{i,8}]
p=Sum[Projection[Cos[x],w[i]],{i,8}]//Together

Plot[{p,Cos[x]},{x,0,Pi}]

Problem 2

Clear[x,v,w]
Do[v[i]=x^(i-1),{i,7}]
iprod[u_,v_]:=Integrate[u*v,{x,-1,1}]
Projection[u_,v_]:=iprod[u,v]/iprod[v,v]v
Do[w[i]=v[i]-Sum[Projection[v[i],w[j]],{j,i-1}],{i,7}]
Do[Print[w[i]//Together],{i,7}]

p=Sum[Projection[Sin[x]+Cos[x],w[i]],{i,7}]//Together

Plot[{p,Sin[x]+Cos[x]},{x,-1,1}]