ME581_Assignment03-Sol

ME581_Assignment03-Sol - ME 581 Assignment #3 - Solution In...

Info iconThis preview shows pages 1–4. Sign up to view the full content.

View Full Document Right Arrow Icon
1 ME 581 Assignment #3 - Solution In the following questions show the derivations that you used to write your program. 1. Write a subroutine to solve the general linear least problem normal equations: } { ] [ } ]{ [ ] [ b A c A A T T = using Cholesky’s [ L ][ L ] T decomposition. The input to your program is matrix [ A ] and a vector { b }. ] [ A : is an n × m matrix ] [ b : is an n -vector ] [ ] [ ] [ A A M T = is a m × m square matrix } { ] [ } { b A d T = is a m- vector We can write the least-squares equations as: } { } ]{ [ d c M = We use Cholesky’s decomposition to find the lower-triangular matrix ] [ L : T L L M ] ][ [ ] [ = Algorithm steps: 1. Calculate: ] [ ] [ ] [ A A M T = 2. Calculate: } { ] [ } { b A d T = 3. Solve } { } ]{ [ d c M = using Cholesky’s T L L ] ][ [ decomposition a. Perform Cholesky’s T L L ] ][ [ decomposition (see Homework 2, Problem 4) b. Solve } { } { ] ][ [ c x L L T = i. Perform back-substitution using a lower triangular matrix } { } ]{ [ c y L = ii. Perform back-substitution using the transpose of a lower triangular matrix } { } { ] [ y x L T = void StandardPartialLeastSquares( double *A, double *b, double *c, int Arows, int Acols) { double M[30*30]; MatrixMultiplyTM(A,Arows,Acols,A,Acols,M); double d[30]; MatrixMultiplyTM(A,Arows,Acols,b,1,d); LL_SolveLinearEquation(M,d,c,Acols); } #define s1(i,j,Mclns) ((i)*(Mclns)+j) //A_Transpose * B void MatrixMultiplyTM( double *A, int Arows, int Acols, double *B, int Bcols, double *C) { int i, j, k, kk; for (i=0;i<Acols;i++){ for (j=0;j<Bcols;j++){ kk=s1(i,j,Bcols); C[kk]=0.; for (k=0;k<Arows;k++)C[kk]+=A[s1(k,i,Acols)]*B[s1(k,j,Bcols)]; } } } void LL_SolveLinearEquation( double *A, double *C, double *X, int n) { double L[30*30], Y[30]; LLcholeskyFactorization(A,L,n); LowerTriangular_BackSubstituion(L,n,C,Y); LowerTriangularTranspose_BackSubstituion(L,n,Y,X); }
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
2 char LLcholeskyFactorization( double *A, double *L, int n) { int i, j, k, nn=n*n; for (i=0;i<nn;i++)L[i]=0.; double Sum=0.; for (i=0;i<n;i++){ for (j=0;j<=i;j++){ Sum=A[s1(i,j,n)]; for (k=0;k<j;k++)Sum=Sum-L[s1(i,k,n)]*L[s1(j,k,n)]; if (i==j)L[s1(i,j,n)]=sqrt(Sum); else L[s1(i,j,n)]=Sum/L[s1(j,j,n)]; } } return (0); } //Back-substitution using a lower triangular matrix void LowerTriangular_BackSubstituion( double *A, int n, double *C, double *X) { int k, j; for (k=0;k<n;k++){ X[k]=0.; for (j=0;j<k;j++)X[k]+=A[s1(k,j,n)]*X[j]; X[k]=(C[k]-X[k])/A[s1(k,k,n)]; } } //Back-substitution using the transpose of a lower triangular matrix void LowerTriangularTranspose_BackSubstituion( double *A, int n, double *C, double *X) { int k, j; for (k=n-1;k>=0;k--){ X[k]=0.; for (j=k+1;j<n;j++)X[k]+=A[s1(j,k,n)]*X[j]; X[k]=(C[k]-X[k])/A[s1(k,k,n)]; } }
Background image of page 2
3 2. Use the subroutine that you wrote in question #1 to find the equation of a parabola that would best fit the equation ) sin( ) ( x x Y = in the interval from x = 0 to x = π . Calculate an estimate of the error between the parabola and the original function using the norm of the residual (error) vector divided by the total number of points.
Background image of page 3

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Image of page 4
This is the end of the preview. Sign up to access the rest of the document.

This note was uploaded on 10/16/2011 for the course MECHANICAL 581 taught by Professor Wasfy during the Fall '11 term at IUPUI.

Page1 / 11

ME581_Assignment03-Sol - ME 581 Assignment #3 - Solution In...

This preview shows document pages 1 - 4. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online