Jeremy Knight Typesetting:-RuleAssistant(); math 639 3/2/11 Algorithm 9.6 - QR Method To obtain the eigenvalues of the symmetric, tridiagonal LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSgmdGltZXM7RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjIyMjIyMjJlbUYnL0ZORlNGK0Y5 matrix LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2J0YrLUkobWFjdGlvbkdGJDYkLUYjNipGKy1GIzYnRistRjI2JC1JKG1mZW5jZWRHRiQ2KS1GIzYnRistSSdtdGFibGVHRiQ2OS1JJG10ckdGJDYqLUkkbXRkR0YkNigtRiM2Ji1JKG1zdWJzdXBHRiQ2KC1GLDYlUSJhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNictSSVtc3VwR0YkNiUtSSNtbkdGJDYkUSIxRicvRlRRJ25vcm1hbEYnLUYjNidGK0ZQLyUrZXhlY3V0YWJsZUdRJmZhbHNlRicvJTBmb250X3N0eWxlX25hbWVHUSgyRH5NYXRoRidGUy8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGUEZdb0Zgb0ZTLUYjNiktSSNtb0dGJDYtUSIoRidGaW4vJSZmZW5jZUdGUi8lKnNlcGFyYXRvckdGX28vJSlzdHJldGNoeUdGUi8lKnN5bW1ldHJpY0dGX28vJShsYXJnZW9wR0Zfby8lLm1vdmFibGVsaW1pdHNHRl9vLyUnYWNjZW50R0Zfby8lJ2xzcGFjZUdRLDAuMTY2NjY2N2VtRicvJSdyc3BhY2VHRlxxRmVuLUZpbzYtUSIpRidGaW5GXHBGXnBGYHBGYnBGZHBGZnBGaHBGanBGXXFGUEZdb0Zgb0ZTRmNvLyUvc3Vic2NyaXB0c2hpZnRHRmVvL0krbXNlbWFudGljc0dGJFEsW25vbmUsbm9uZV1GJ0Zdby8lLHBsYWNlaG9sZGVyR0ZSRmluLyUpcm93YWxpZ25HRi4vJSxjb2x1bW5hbGlnbkdGLi8lK2dyb3VwYWxpZ25HRi4vJShyb3dzcGFuR1EiMUYnLyUrY29sdW1uc3BhbkdGYXItRkY2KC1GIzYnLUZLNigtRiw2JVEiYkYnRlBGUy1GIzYnLUZmbjYkUSIyRidGaW5GUEZdb0Zgb0ZTRmZvRmNvRmJxRmRxRitGXW9GZ3FGaW5GaXFGW3JGXXJGX3JGYnItRkY2KC1GIzYmLUZmbjYkUSIwRidGaW5GXW9GZ3FGaW5GaXFGW3JGXXJGX3JGYnItRkY2KC1GIzYoRistRiM2JS1GLDYlUSMuLkYnRlBGU0Zdb0Zpbi1GaW82LVEiLkYnRmluL0ZdcEZfb0ZecC9GYXBGX29GYnBGZHBGZnBGaHAvRltxUSYwLjBlbUYnL0ZecUZodEZdb0ZncUZpbkZpcUZbckZdckZfckZickZic0ZpcUZbckZdci1GQzYqLUZGNigtRiM2JkZockZdb0ZncUZpbkZpcUZbckZdckZfckZici1GRjYoLUYjNiYtRks2KEZNRl1zRmZvRmNvRmJxRmRxRl1vRmdxRmluRmlxRltyRl1yRl9yRmJyLUZGNigtRiM2Ji1GaW82LVEoJmR0ZG90O0YnRmluRmV0Rl5wRmZ0RmJwRmRwRmZwRmhwRmd0Rml0Rl1vRmdxRmluRmlxRltyRl1yRl9yRmJyRmZ1LUZGNigtRiM2Ji1GaW82LVEpJnZlbGxpcDtGJ0ZpbkZldEZecEZmdEZicEZkcEZmcEZocEZndEZpdEZdb0ZncUZpbkZpcUZbckZdckZfckZickZpcUZbckZdci1GQzYqRmJzRmZ1RmZ1RmZ1RmJzRmlxRltyRl1yLUZDNipGXXZGZnVGZnVGZnUtRkY2KC1GIzYmLUZLNihGanItRiM2Jy1GLDYlUSJuRidGUEZTRlBGXW9GYG9GU0Zmb0Zjb0ZicUZkcUZdb0ZncUZpbkZpcUZbckZdckZfckZickZpcUZbckZdci1GQzYqRmJzLUZGNigtRiM2Ji1GaW82LVEoJmN0ZG90O0YnRmluRmV0Rl5wRmZ0RmJwRmRwRmZwRmhwRmd0Rml0Rl1vRmdxRmluRmlxRltyRl1yRl9yRmJyRmJzRmh2LUZGNigtRiM2Ji1GSzYoRk1GXndGZm9GY29GYnFGZHFGXW9GZ3FGaW5GaXFGW3JGXXJGX3JGYnJGaXFGW3JGXXIvJSZhbGlnbkdRJWF4aXNGJy9GanFRKWJhc2VsaW5lRicvRlxyUSdjZW50ZXJGJy9GXnJRJ3xmcmxlZnR8aHJGJy8lL2FsaWdubWVudHNjb3BlR0ZSLyUsY29sdW1ud2lkdGhHUSVhdXRvRicvJSZ3aWR0aEdGX3kvJStyb3dzcGFjaW5nR1EmMS4wZXhGJy8lLmNvbHVtbnNwYWNpbmdHUSYwLjhlbUYnLyUpcm93bGluZXNHUSVub25lRicvJSxjb2x1bW5saW5lc0dGankvJSZmcmFtZUdGankvJS1mcmFtZXNwYWNpbmdHUSwwLjRlbX4wLjVleEYnLyUqZXF1YWxyb3dzR0Zfby8lLWVxdWFsY29sdW1uc0dGX28vJS1kaXNwbGF5c3R5bGVHRl9vLyUlc2lkZUdRJnJpZ2h0RicvJTBtaW5sYWJlbHNwYWNpbmdHRmd5RitGXW9GaW4vJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJ0Zpbi9GZXFRJ01hdHJpeEYnLyUlb3BlbkdRIltGJy8lJmNsb3NlR1EiXUYnRmBbbC8lK2FjdGlvbnR5cGVHUS5ydGFibGVhZGRyZXNzRidGK0Zdb0ZpbkYrRlBGXVtsRl1vRmBvRlNGaFtsRitGXW9GaW5GK0Zdb0Zpbg==
<Text-field style="Heading 1" layout="Heading 1">QR(n,a0,b0,TOL,M)</Text-field> QR:=proc(n0,a0,b0,TOL,M) #n = dimension; a = 1xn array from diagonal; #b= 1xn array with 0 as the first entry from the 2nd diagonal #TOL= tolerance; M = maximum number of iterations. local a,b,c,d,b1,c1,d1,i,j,k,n,x,y,z,q,r,mu1,mu2,sigma,sigma1,shift, lambda,lambda1,lambda2,a_split1,a_split2,b_split1,b_split2,temp1,temp2,eigenvalues: #option trace: #Check input to require b0 to have 0 as first entry if b0[1]<>0 then print ("b needs to be a 1xn array with 0 as the first entry"); break: end if: #Initialize variables a:=a0: b:=b0: n:=n0: k:=1: shift:=0: #Initialize eigenvalues array to store the eigenvalues eigenvalues:=Array(1..n); ### Main loop while k<=M do if abs(b[n]) <= TOL then lambda:=a[n]+shift: print("lambda",n, lambda): eigenvalues[n]:=lambda: n:=n-1: end if: if abs(b[2])<=TOL then lambda:=a[1]+shift: print("lambda",n, lambda): eigenvalues[n]:=lambda: n:=n-1: a[1]:=a[2]: for j from 2 to n do a[j]:=a[j+1]: b[j]:=b[j+1]: end do: end if: if n=0 then break: end if: if n=1 then lambda:=a[1]+shift: print("lambda",n, lambda): eigenvalues[n]:=lambda: break: end if: #if b[j]<=TOL we will define the diagonals of the split matrices for j from 3 to n-1 do if b[j]<=TOL then a_split1:=Array(1..(j-1)): b_split1:=Array(1..(j-1)): a_split2:=Array(1..(n-j+1)): b_split2:=Array(1..(n-j+1)): for i from 1 to j-1 do a_split1[i]:=a[i]: b_split1[i]:=b[i]: end do: for i from j to n do a_split2[i]:=a[i]: if i=j then b_split2[i]:=0: else b_split2[i]:=b[i]: end if: end do: print("split into the following vectors a_i,b_i"): print(a_split1):print(b_split1): print("and"): print(a_split1):print(b_split2): end if: end do: # Compute Shift b1:=-(a[n-1]+a[n]): c1:=a[n]*a[n-1]-(b[n])^2: d1:=(b1^2-4.0*c1)^(1/2): if b1 > 0 then mu1:=(-2*c1)/(b1+d1): mu2:=-(b1+d1)/2: else mu1:=(d1-b1)/2: mu2:=(2*c1)/(d1-b1): end if: if n=2 then lambda1:=mu1+shift: lambda2:=mu2+shift: print("lambda[1] and lambda[2] =",lambda1,lambda2): eigenvalues[1]:=lambda1: eigenvalues[2]:=lambda2: break: end if: #compute sigma1 to use in shift temp1:=abs(mu1-a[n]): temp2:=abs(mu2-a[n]): if temp1<temp2 then sigma1:=mu1: else sigma1:=mu2: end if: #Accumulate the shift shift:=shift+sigma1: #Perform Shift d:=Array(1..n): for j from 1 to n do d[j]:=a[j]-sigma1: end do: #Compute R^(k) and A^(k)=P_j*A_(j-1) x:=Array(1..n): y:=Array(1..n): sigma:=Array(1..n): x[1]:=d[1]: y[1]:=b[2]: for j from 2 to n do z[j-1]:=((x[j-1])^2 + (b[j])^2)^(1/2): c[j]:=x[j-1]/z[j-1]: sigma[j]:=b[j]/z[j-1]: q[j-1]:=c[j]*y[j-1]+sigma[j]*d[j]: x[j]:=-sigma[j]*y[j-1]+c[j]*d[j]: if j<>n then r[j-1]:=sigma[j]*b[j+1]: y[j]:=c[j]*b[j+1]: end if: end do: #Compute A^(k+1) z[n]:=x[n]: a[1]:=sigma[2]*q[1]+c[2]*z[1]: b[2]:=sigma[2]*z[2]: for j from 2 to n-1 do a[j]:=sigma[j+1]*q[j]+c[j]*c[j+1]*z[j]: b[j+1]:=sigma[j+1]*z[j+1]: end do: a[n]:=c[n]*z[n]: k:=k+1: a; b; end do: ###end main loop if k>M then print("Maximum number of iterations exceeded"): end if: return eigenvalues; end proc:
JSFH LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= a:=Array([3.0,3.0,3.0]); b:=Array([0.0,1.0,1.0]); QR(3,a,b,.000001,10); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIzo9RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRj0vJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRkwtSShtYWN0aW9uR0YkNiUtSShtZmVuY2VkR0YkNigtSSdtdGFibGVHRiQ2NS1JJG10ckdGJDYoLUkkbXRkR0YkNigtSSNtbkdGJDYkUSQzLjBGJ0Y5LyUpcm93YWxpZ25HUSFGJy8lLGNvbHVtbmFsaWduR0Zeby8lK2dyb3VwYWxpZ25HRl5vLyUocm93c3BhbkdRIjFGJy8lK2NvbHVtbnNwYW5HRmVvRmVuRmVuRlxvRl9vRmFvLyUmYWxpZ25HUSVheGlzRicvRl1vUSliYXNlbGluZUYnL0Zgb1EnY2VudGVyRicvRmJvUSd8ZnJsZWZ0fGhyRicvJS9hbGlnbm1lbnRzY29wZUdGMS8lLGNvbHVtbndpZHRoR1ElYXV0b0YnLyUmd2lkdGhHRmVwLyUrcm93c3BhY2luZ0dRJjEuMGV4RicvJS5jb2x1bW5zcGFjaW5nR1EmMC44ZW1GJy8lKXJvd2xpbmVzR1Elbm9uZUYnLyUsY29sdW1ubGluZXNHRmBxLyUmZnJhbWVHRmBxLyUtZnJhbWVzcGFjaW5nR1EsMC40ZW1+MC41ZXhGJy8lKmVxdWFscm93c0dGPS8lLWVxdWFsY29sdW1uc0dGPS8lLWRpc3BsYXlzdHlsZUdGPS8lJXNpZGVHUSZyaWdodEYnLyUwbWlubGFiZWxzcGFjaW5nR0ZdcUY5L0krbXNlbWFudGljc0dGJFEnVmVjdG9yRicvJSVvcGVuR1EiW0YnLyUmY2xvc2VHUSJdRidGY3IvJSthY3Rpb250eXBlR1EucnRhYmxlYWRkcmVzc0YnLyUpcnRhYmxlaWRHUSoyNDk5NjY4MzJGJy8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJSlyZWFkb25seUdGPUY5 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEiYkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIzo9RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRj0vJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRkwtSShtYWN0aW9uR0YkNiUtSShtZmVuY2VkR0YkNigtSSdtdGFibGVHRiQ2NS1JJG10ckdGJDYoLUkkbXRkR0YkNigtSSNtbkdGJDYkUSMwLkYnRjkvJSlyb3dhbGlnbkdRIUYnLyUsY29sdW1uYWxpZ25HRl5vLyUrZ3JvdXBhbGlnbkdGXm8vJShyb3dzcGFuR1EiMUYnLyUrY29sdW1uc3BhbkdGZW8tRmZuNigtRmluNiRRJDEuMEYnRjlGXG9GX29GYW9GY29GZm9GaG9GXG9GX29GYW8vJSZhbGlnbkdRJWF4aXNGJy9GXW9RKWJhc2VsaW5lRicvRmBvUSdjZW50ZXJGJy9GYm9RJ3xmcmxlZnR8aHJGJy8lL2FsaWdubWVudHNjb3BlR0YxLyUsY29sdW1ud2lkdGhHUSVhdXRvRicvJSZ3aWR0aEdGanAvJStyb3dzcGFjaW5nR1EmMS4wZXhGJy8lLmNvbHVtbnNwYWNpbmdHUSYwLjhlbUYnLyUpcm93bGluZXNHUSVub25lRicvJSxjb2x1bW5saW5lc0dGZXEvJSZmcmFtZUdGZXEvJS1mcmFtZXNwYWNpbmdHUSwwLjRlbX4wLjVleEYnLyUqZXF1YWxyb3dzR0Y9LyUtZXF1YWxjb2x1bW5zR0Y9LyUtZGlzcGxheXN0eWxlR0Y9LyUlc2lkZUdRJnJpZ2h0RicvJTBtaW5sYWJlbHNwYWNpbmdHRmJxRjkvSSttc2VtYW50aWNzR0YkUSdWZWN0b3JGJy8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJ0Zoci8lK2FjdGlvbnR5cGVHUS5ydGFibGVhZGRyZXNzRicvJSlydGFibGVpZEdRKjI1MDI5OTM2NEYnLyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GJy8lKXJlYWRvbmx5R0Y9Rjk= LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbXNHRiQ2I1EnbGFtYmRhRictSSNtb0dGJDYtUSIsRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR1EldHJ1ZUYnLyUpc3RyZXRjaHlHRjgvJSpzeW1tZXRyaWNHRjgvJShsYXJnZW9wR0Y4LyUubW92YWJsZWxpbWl0c0dGOC8lJ2FjY2VudEdGOC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUkjbW5HRiQ2JFEiM0YnRjNGLy1GTTYkUSwxLjU4NTc4NjQzOEYnRjNGMw== LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbXNHRiQ2I1E6bGFtYmRhWzFdfmFuZH5sYW1iZGFbMl1+PUYnLUkjbW9HRiQ2LVEiLEYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdRJXRydWVGJy8lKXN0cmV0Y2h5R0Y4LyUqc3ltbWV0cmljR0Y4LyUobGFyZ2VvcEdGOC8lLm1vdmFibGVsaW1pdHNHRjgvJSdhY2NlbnRHRjgvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1JI21uR0YkNiRRLDQuNDE0MjEzNTYyRidGM0YvLUZNNiRRLDMuMDAwMDAwMDAwRidGM0Yz LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKm9JSV0j LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= TTdSMApJNlJUQUJMRV9TQVZFLzI0OTk2NjgzMlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKzpEcip6IyEiKiQiK2VfKUdXIkYpJCErSHN0WEAhIzk2Ig==TTdSMApJNlJUQUJMRV9TQVZFLzI1MDI5OTM2NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiIiFGKCQiKzdWJ1EqPiEjNSQhK0hId3pXISM7NiI=TTdSMApJNlJUQUJMRV9TQVZFLzI1MDMwMzA2OFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK2lOQDlXISIqJCIrKysrK0lGKSQiK1FreSZlIkYpNiI=