{VERSION 2 3 "IBM INTEL NT" "2.3" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 12 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }} {SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 0 "" }}{PARA 18 "" 0 "" {TEXT -1 16 "The buckled beam" }}{PARA 19 "" 0 "" {TEXT -1 12 "by Mary Pugh" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 92 " As described in class, we are trying to find how much force can be app lied to beam until it " }}{PARA 0 "" 0 "" {TEXT -1 96 "buckles. The b eam is of length L, with x = 0 being its left end, and x = L being its right end." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 95 "The height of the beam is given by y(x), with the flat beam bei ng y(x) = 0 for all x. The ends" }}{PARA 0 "" 0 "" {TEXT -1 43 "of th e beams are fixed: y(0) = 0, y(L) = 0." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 54 "The beam satisfies the equation \+ R y_xx = - P y" }}{PARA 0 "" 0 "" {TEXT -1 70 "where R is the flexu ral rigidity and P is the load placed on the beam." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 91 "You should imagine a beam sticking horizontally out of the wall. You push on the beam and " }} {PARA 0 "" 0 "" {TEXT -1 89 "if you push weakly, the beam stays straig ht. As you push harder, the beam should buckle." }}{PARA 0 "" 0 "" {TEXT -1 87 "At some point, the beam will snap. This involves complic ated nonlinear effects, which " }}{PARA 0 "" 0 "" {TEXT -1 31 "the equ ation does not consider." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 99 "The question is: for what loads are there nontrivi al solutions? If the load is weak, then the only" }}{PARA 0 "" 0 "" {TEXT -1 95 "solution should be the straight beam y(x) = 0. If the lo ad reaches a critical value, the beam " }}{PARA 0 "" 0 "" {TEXT -1 93 "should buckle and y(x) will be some function other than 0. If the lo ad gets even larger, the" }}{PARA 0 "" 0 "" {TEXT -1 31 "beam should b uckle yet again. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 93 "Okay, so we have some differential equation we want to so lve. We will approximate a solution" }}{PARA 0 "" 0 "" {TEXT -1 101 " using \"finite difference methods\". (In math 320, you learn more abo ut putting differential equations" }}{PARA 0 "" 0 "" {TEXT -1 15 "on a computer.)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 50 "We consider N+1 equally spaced points on the beam:" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 103 " \+ 0 = x_0 < x_1 < x_2 < x_3 < ...... < x_(n-2) < x _(N-1) < x_N = L" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 53 "So we'll ask that the equation holds at these points:" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 103 " \+ R y_xx(x_i) = - P y(x_i) \+ for i=0,1,...N" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 101 "Now. If you only have the solution at N+1 points , how can you use these N+1 pieces of information to" }}{PARA 0 "" 0 " " {TEXT -1 39 "approximate the second derivative of y?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 86 "You know that dy/d x = lim (y(x+h)-y(x))/h as h goes to zero. So we use this to say:" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 97 " \+ dy/dx (x_i) ~ ( y(x_(i+1) ) - y(x_i) )/delta" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 67 "where delta is the distance between points: delta = x_ (i+1)-x_i." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 "We do this again, and approximate the second derivative:" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 81 " \+ d^2y/dx^2 = d/dx dy/dx ~ ( dy/dx(x_i) - dy/dx(x_(i-1) )/delt a " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 104 "We plug in dy/dx (x_i) ~ ( y(x_(i+1)) - y(x_i) )/delta dy/dx (x_i) \+ ~ ( y(x_i) - y(x_(i-1)) )/delta" }}{PARA 0 "" 0 "" {TEXT -1 11 "This \+ yields" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 80 " d^2y/dx^2 (x_i) ~ (y_(i+1) - 2 y_i + y_(i- 1))/delta^2" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 66 "So the equation R y_xx = - P y is approximated with N-1 equatio ns:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 61 " \+ R (y_2 - 2 y_1)/delta^2 = - P y_1" }}{PARA 0 "" 0 "" {TEXT -1 59 " R(y_3 - 2 y_2 + y_1)/delta ^2 = - P y_2" }}{PARA 0 "" 0 "" {TEXT -1 59 " R(y_ 4 - 2 y_3 + y_2)/delta^2 = - P y_3" }}{PARA 0 "" 0 "" {TEXT -1 47 " \+ ........" }}{PARA 0 "" 0 "" {TEXT -1 65 " R(y_(n-1) - 2 y_(n-2) + y_(n-3))/delta^2 = - P y_(n- 2)" }}{PARA 0 "" 0 "" {TEXT -1 65 " R( - 2 y_(n-1) + \+ y_(n-2))/delta^2 = - P y_(n-1)" }}{PARA 0 "" 0 "" {TEXT -1 1 " " }} {PARA 0 "" 0 "" {TEXT -1 108 "In the above, I'm using short-hand: y_ i = y(x_i). Also, I have used the fact that y_0 = y(x_0) = y(0) = 0" }}{PARA 0 "" 0 "" {TEXT -1 28 "and y_N = y(x_N) = y(L) = 0." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 99 "So we have N-1 \+ linear equations, involving N-1 unknowns y_1, ... , y_(N-1). We write the equations" }}{PARA 0 "" 0 "" {TEXT -1 68 "as A y = - P y wher e A is a matrix with the following structure:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 105 " A is -2 R/delta^ 2 on the diagonal, and R/delta^2 right above and below the diagonal, 0 elsewhere." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 81 "We are trying to find the eigenvalues for the problem. So to d o this we need to " }}{PARA 0 "" 0 "" {TEXT -1 31 "1) choose a length \+ for the beam" }}{PARA 0 "" 0 "" {TEXT -1 90 "2) choose the number of p oints --- this determines both delta and the size of the matrix A" }} {PARA 0 "" 0 "" {TEXT -1 45 "3) choose a value for the flexural rigidi ty R" }}{PARA 0 "" 0 "" {TEXT -1 22 "4) define the matrix A" }}{PARA 0 "" 0 "" {TEXT -1 26 "5) compute its eigenvalues" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(linalg ):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 91 "We first pick the number of points and the length of the beam. This then determines delta," }} {PARA 0 "" 0 "" {TEXT -1 67 "the distance between points. We also cho ose the flexural rigidity." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "N := \+ 10; L := 3; delta := L/N; R := 2;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#> %\"NG\"#5" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"LG\"\"$" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#>%&deltaG#\"\"$\"#5" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"RG\"\"#" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "We \+ now define the matrix A:" }{MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "A := matrix(N-1,N-1): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to N-1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 " \+ for j from 1 to N-1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 " A[i, j] := 0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 " od" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to N-1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 " A[i,i] := -2*R/d elta^2:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to N-2 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 " A[i,i+1] := R/delta^2:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 " A[i+1,i] := R/delta^2:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "o d:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "print(A): " }}{PARA 0 "" 0 " " {TEXT -1 81 "Note: if you don't want to see the matrix, remove the \+ \"print(A)\" statement above." }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'MA TRIXG6#7+7+#!$+%\"\"*#\"$+#F*\"\"!F-F-F-F-F-F-7+F+F(F+F-F-F-F-F-F-7+F- F+F(F+F-F-F-F-F-7+F-F-F+F(F+F-F-F-F-7+F-F-F-F+F(F+F-F-F-7+F-F-F-F-F+F( F+F-F-7+F-F-F-F-F-F+F(F+F-7+F-F-F-F-F-F-F+F(F+7+F-F-F-F-F-F-F-F+F(" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "Now we compute the eigenvalues:" }{MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "lambda := eval f(Eigenvals(A));" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%'lambdaG-%'VECTO RG6#7+$!+%Hi8n)!\")$!+Hb2S!)F+$!+UL#o0(F+$!+1L&y\"eF+$!+XWWWWF+$!+#eN5 2$F+$!+Vb1K=F+$!+!eL\")[)!\"*$!+_fEv@F:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 87 "The matrix A is symmetric and real, hence by the theorem, all the eigenvalues are real." }} {PARA 0 "" 0 "" {TEXT -1 95 "(Also, since the matrix is symmetric, it \+ is easier to find the eigenvalues with a computer than" }}{PARA 0 "" 0 "" {TEXT -1 33 "if the matrix weren't symmetric.)" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 84 "Thus we have found the \+ eigenvalues. These will satisfy A v = lambda v. Since we're" }} {PARA 0 "" 0 "" {TEXT -1 85 "looking for solutions of A y = - P y, we \+ are interested in the negatives of the above" }}{PARA 0 "" 0 "" {TEXT -1 5 "list." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 67 "Okay, let's think. We can acutally solve the problem analytica lly:" }}{PARA 0 "" 0 "" {TEXT -1 61 " R y_xx = - P y with y(0) = 0, y(L) = 0" }}{PARA 0 "" 0 "" {TEXT -1 84 "suggests \+ a solultion of y(x) = A sin(k x Pi/L) since y(0) = y(L) = 0. Plugging this" }}{PARA 0 "" 0 "" {TEXT -1 24 "into the equation yields" }} {PARA 0 "" 0 "" {TEXT -1 69 " - k^2 pi^2/L^2 R A sin(k x Pi/L) = - P A sin(k x Pi/L)" }}{PARA 0 "" 0 "" {TEXT -1 87 "this will always have a zero solution: A = 0. But to have a nonzero solution, \+ we would" }}{PARA 0 "" 0 "" {TEXT -1 12 "have to have" }}{PARA 0 "" 0 "" {TEXT -1 64 " - k^2 Pi^2/L^2 R = lamd a = - P " }}{PARA 0 "" 0 "" {TEXT -1 30 "This is something we can te st:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "lambda[N-1];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "k \+ := 1; test := evalf(-k^2*Pi^2*R/L^2);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "lambda[N-1]-test;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$!+_fEv@!\" *" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"kG\"\"\"" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%testG$!+BaC$>#!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 " 6#$\")r%zz\"!\"*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "lambda[ N-2];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "k := 2; test := evalf(-k^2 *Pi^2*R/L^2);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "lambda[N-2]-test; " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$!+!eL\")[)!\"*" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#>%\"kG\"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%tes tG$!+$p\")Hx)!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"*8\"[[G!\"*" } }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 93 " Note: we're lucky this time around that the equation is exactly solvab le. Symmetric matrices" }}{PARA 0 "" 0 "" {TEXT -1 92 "show up all ov er the place, and very often one can't find exact solutions. This mea ns that " }}{PARA 0 "" 0 "" {TEXT -1 84 "either you put the problem on a computer or you try to do something clever involving" }}{PARA 0 "" 0 "" {TEXT -1 95 "approximating the solutions. In both cases, you'll \+ be _very_ thankful the matrix is symmetric." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 11 "Problem 1) " }}{PARA 0 "" 0 "" {TEXT -1 88 " Take your social security number. Add up its digits. D ivide by 10. Take this as the " }}{PARA 0 "" 0 "" {TEXT -1 22 " lengt h of the bar, L." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 95 " Graph the smallest eigenvalue as a function of the flexu ral rigidity. The smallest eigenvalue" }}{PARA 0 "" 0 "" {TEXT -1 90 " determines the largest load that can be put on the beam before it bu ckles. How does this" }{MPLTEXT 1 0 0 "" }}{PARA 0 "" 0 "" {TEXT -1 84 " critical load depend on the rigidity? Explain whether this makes sense physically." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 10 "Problem 2)" }}{PARA 0 "" 0 "" {TEXT -1 89 " Take your so cial security number. Add up its digits. Divide by 10. Take this as the " }}{PARA 0 "" 0 "" {TEXT -1 56 " length of the bar, L. Take fl exural rigidity R = 1. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 70 " Now use the above routine to compute the critica l load as follows: " }}{PARA 0 "" 0 "" {TEXT -1 78 " go to the top o f the worksheet. Modify L and R appropriately. Set N to 10." }} {PARA 0 "" 0 "" {TEXT -1 87 " Execute the commands, building the matr ix, finding its eigenvalues, and comparing the" }}{PARA 0 "" 0 "" {TEXT -1 84 " computed eigenvalues to the true answers. (Only do the two smallest eigenvalues.)" }}{PARA 0 "" 0 "" {TEXT -1 83 " Make a t able which has: N, error on first eigenvalue, error on second eigenval ue." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 " \+ Go back to the top of the worksheet. Set N to 20. " }}{PARA 0 "" 0 " " {TEXT -1 87 " Execute the commands, building the matrix, finding it s eigenvalues, and comparing the" }}{PARA 0 "" 0 "" {TEXT -1 84 " com puted eigenvalues to the true answers. (Only do the two smallest eige nvalues.)" }}{PARA 0 "" 0 "" {TEXT -1 78 " Add to your table: N, erro r on first eigenvalue, error on second eigenvalue." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 " Go back to the top of t he worksheet. Set N to 30. " }}{PARA 0 "" 0 "" {TEXT -1 87 " Execut e the commands, building the matrix, finding its eigenvalues, and comp aring the" }}{PARA 0 "" 0 "" {TEXT -1 84 " computed eigenvalues to th e true answers. (Only do the two smallest eigenvalues.)" }}{PARA 0 " " 0 "" {TEXT -1 78 " Add to your table: N, error on first eigenvalue, error on second eigenvalue." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 " Go back to the top of the worksheet. Set N t o 40. " }}{PARA 0 "" 0 "" {TEXT -1 87 " Execute the commands, buildi ng the matrix, finding its eigenvalues, and comparing the" }}{PARA 0 " " 0 "" {TEXT -1 84 " computed eigenvalues to the true answers. (Only do the two smallest eigenvalues.)" }}{PARA 0 "" 0 "" {TEXT -1 78 " A dd to your table: N, error on first eigenvalue, error on second eigenv alue." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 95 " Keep increasing N by 10 until it you lose patience waiting for it to compute the eigenvalues." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 85 " How are the errors behaving as N increases? Can \+ you make a guess at a function of N" }}{PARA 0 "" 0 "" {TEXT -1 53 " t hat might describe the errors? I.e., error = F(N)." }}}}{MARK "8 4 0 " 95 }{VIEWOPTS 1 1 0 1 1 1803 }