MIT OpenCourseWare http://ocw.mit.edu 18.085 Computational Science and Engineering IFall 2008 For information about citing these materials or our Terms of Use, visit: http://ocw.mit.edu/terms.18.085: Matlab Homework #5 Laplace's equation. The problem is Laplace's eqixa,tio~i on the unit sqiiare with boiinda,ry condi- tions u= y on the side :r = 1?u= n: on the side y = 1, 117 .n, = -y on t,he side n: =0, 111 : n,= -n: on the side y = 0. Replace the 2nd derivatives in Laplace's equation by centered second differ- ences. This gives the "5-point discrete L.aplacianV on a sqixsre grid (An: = Ay): [u(:I:+AT, g) +IL(Z-Ax, y) +U(X,?/+Ax) +u (x: -Ax) -44r, ?I)]/ (4~)~ Set Aa = 1/11 giving 10 x 10 = 100 int,erior grid points with 100 unknowns ~(r, The 100 grid points fall in a sqiia,re array hiit yon have t,o make 11). them into a VECTOR with 100 components. I usually number them by romTs;11.1 = AX, Ax) and uz = 11,(Ax, 2Az) and next row 11,11 = 11,(2An:, Ax) and last corner uloo = u(1OAx. 10Az). At turo boundaries ure know u.At t,he left boixnda,ry n: = 0 ure know II! = -?I. Replace by [~r(Ax, y) -u(0, ?/)]/Ax = -y. This gives lr(0, ?I)in terms of u(Az, Ax); Siibstitut,e in the 5-point eqiiation t,o elimi~iate u(0,y) . Simi-larly, eliminate u,(:r, 0) on the bolxnda,ry y = 0 where do,/&/ = -x. Set up t,he whole system as KT],= f where K is 100 x 100. Is K symmet-ric? Is positive definite?? (Let k4ATL.AB decide.) Print out diag(K). IT SH0UL.D NOT BE ALL 4's. K is ba~ided armind the ma,in diagonal. What, is t,he bandwidt,h so K(i, j) = 0 if abs (i-j) > IN? FVha,t are the largest a~id smallest eigenvahxes of K? Cse the command eig(K). Print out 11, the first coliimn or row of K inverse. This gives the valixe of AX, Ax) at the lower left corner from the right side f. What is the ratio of the first compo~ient v(l? 1) to the last component ~(10, lo)? Bonus: Let Matlab solve Ku = 0 a,nd plot u(n:,;v) (Remember, 71,(z,?j) is a 2D sixrface). Mattlab also has a plotting command called QUIVER which plots velocity fields. If u(a,y) describes pot,ent,ial flow, there is a ve1ocit;y field associated with ?I, given by v = Vu. Plot this velocit,y
View Full Document
Unlocking...