HW5 Solution

5.1 Fundamentals: Jacobi and Gauss-Seidel Methods

3 pts Consider the 5 equa;ons for 5 unknowns, wriKen in matrix form at right.

Reorder the equa;ons to form a new A x = b problem, wriKen in matrix form, where the new matrix A is “strictly diagonally dominant” (or at least the “best you can do” to make as “strong” a diagonal as possible).



−2
0
3
−6
−1

5
6
7
8
9







2
4
−4
1
2

7
4
−1
7
2







−1
3
−1
0
0






x
1
x
2
x
3

x4
x
5








5







4






=
3







2











1























5.2 Fundamentals: Jacobi and Gauss-Seidel Methods

5 pts Consider the 4x4 problem A x = b at right, for which I have already ordered equa;ons to create as “strong” a diagonal as I can (at least the largest element in each row is on the diagonal).















x





2
−1 0
1


1



1





x






1
−2 1
1


2



1








=


−114


x


−1

−1





1
0−12

x
3



2





4
















Star;ng with the ini;al guess vector x0 = [0 1 0 -1]T , calculate by-hand (and show your work!!) ...

+3 a) the first two itera;ons of the solu;on (x1 and x2) using the Jacobi itera;on method.

+2 b) the first two itera;ons of the solu;on (x1 and x2) using the Gauss-Seidel itera;on method.











5.3
Fundamentals: Errors and Norms.







0
















2





6 pts
Consider the same 4x4 problem from 5.2, with the exact solu;on x and


x
=








itera;ons x5 and x6 given at right. Do the following by-hand:







1































−1




+3
a) For the k=6th itera;on, evaluate the true error, the proxy for the

















0.2







0.5


error, and the residual error, expressed as vectors (
don’t
take norms).



































5 =

1.7




6 =

2.3

+3
b) Evaluate the L1, L2, and L∞ norms of the exact solu;on vector x.
x






x





Show your work and explain everything you do



0.7







0.7






0.2







−0.5



















HW5 February 15

PLEASE READ THIS! It sets up the “Sagging Cable” problem 5.4
A flexible cable is supported between two poles, height YL at le]
and YR at right, separated by distance L.
d
2
y

1 dy 1
A differen;al equa;on represen;ng the

















dx2
− 2a dx = a
cable shape can be approximated as:


where a is a constant represen;ng the proper;es of the cable (light, high-tension cables have a high value of a; heavy, low-

tension cables have a low value of a). Given the boundary condi;ons in the figure shown, across poles separated by L = 10 m, and for a cable with property a = 2, you could

YL=20
20





15

y(x)


YR=15
10





5

L=10









0
2
4
6
8
10
0





x=0




x=L

solve this equa;on using calculus from your other math classes, but here we’ll solve it NUMERICALLY.

Our way to solve the problem is to define variables Y to Y
to
Y1










1
M
20





represent the height every x across the space x = 0 to L. Then

Y2



Y6
it turns out you can approximate the deriva;ves above as:
15

Y3

Y5





Y4


d
2
y(xi )
≈ Yi−1 − 2Yi +Yi+1
dy(xi )
≈ −Yi−1 +Yi+1












10






dx2
( x)2
dx
2 x


x
x
x
x
x
(You’ll learn more about approxima;ng deriva;ves in a month
5












or so; for now, just trust me that this approxima;on works!)






Then subs;tu;ng all this into the differen;al equa;on gives:
0
2
4
6
8
10







0






1

1



−2



1

1



1



+

Y
+

Y
+

−

Y
=



( x)2




( x)2


( x)2




a




4a x
i−1



i



4a x
i+1




This linear equa;on applies for all the “interior” points i (i.e. 2 ≤ i ≤ M-1). Coupled with the two boundary
condi;ons segng the le] and right heights (Y1 = YL = 20, and YM = YR = 15) I now can write M equa;ons for
all M unknown Yi !!







Y1=20


ß LeG boundary condn

For Example: Take M = 6 unknown points (Y1 to Y6) across the span,


5


1


3

1


spaced by
x = 2 meters, as shown by the blue dots above. Then,












(i = 2:5)
given cable a = 2, I can write six equa;ons for the unknowns as:

16 Yi−1 −
2 Yi
+ 16 Yi+1 =
2









Finally, I can collect all M = 6 equa;ons in matrix form as:



Y6 =15


ß Right boundary condn


Y












1
0
0
0
0
0




20

ß LeG boundary cond
n








1






516 −12 316
0
0
0

Y

1 2














2


































0
5 16
−1 2
3 16
0
0

Y3


1 2



Interior equaLons

0
0
5 16
−1 2
3 16
0

Y

=
1 2



for 2 ≤ i ≤ 5






























4





















Y5














0
0
0
516 −12 316




1 2










0
0
0
0
0
1




15

















Y




ß Right boundary condn






























6














No;ce that I ordered the equa;ons the way I did to make the diagonal as strong as possible!

Your job in problem 5.4 is to extend this same example to many more points, by crea;ng the matrices for M = 201 equa;ons for unknown variables Yi to Y201, and then solving for all the Yi using Jacobi itera;on.

HW5











February 15
5.4
Jacobi IteraTon MATLAB code




Y1




12 pts
Extend the concept from the last page to use the Jacobi
20








Y3



YM


Itera;on method to solve for the height Y1 through YM
15









Y2

…
YM-2



at M = 201 equally-spaced (
x) points along a cable.





YM-1










Assume characteris;cs of the cable are the same as before: L = 10 10

x









meters, a = 2, and boundary condi;ons YL = 20 and YR = 15.






Like last page, I can approximate the differen;al equa;on with the 5





linear rela;onship below:








0






1
2 +
1


−2



1
2 −
1

1 0
2
x
6
8
10

( x)


Yi−1
+
( x)
2
Yi +
( x)


Yi+1 =
a




L=10



4a x







4a x







If you don’t really understand exactly where this comes from, don’t worry! You’re going to be deriving this yourself by Chapter 10. For now, just accept the boxed equa;on applies for any three adjacent Yi-1, Yi, Yi+1 (2 ≤ i ≤ M-1) separated by a constant x. The only difference to the example on the last page is that M is now 201 instead of 6, so the x is much smaller (10/200 = 0.05, instead of 2 meters). So, including the two boundary condi;ons (Y1 = YL and YM = YR), you now can write 201 equa;ons for all 201 unknown Yi !!

To solve the problem I have posted a script HW5_4.m that already does all the following for you:

• 

• 

• 

• 

• 

Ini;alizes variables for the rod (L, a, YL, YR), and M = 201,
sets values for convergence tolerances for the L2 norm of “proxy” error between itera;ons and the residual error (tolX = tolR = sqrt(M) x 10-8, using M = 201),
calls a func;on Jacobi.m to itera;vely solve for all the Yi un;l convergence,
^
calculates the exact analy;c solu;on Y( x) (to compare with your itera;ve one),

makes two plots: (i) the Yi(xi) distribu;on for the 1st, 100th, 1000th, 10 000th and final (converged) itera;ons compared to the exact solu;on, and (ii) the log10(error) convergence history.

Your job is just to make the func;on Jacobi.m that does the following:

a) EXACTLY interfaces with my main script, using all the same inputs and outputs. To help, I’ve actually provided the first line of the func;on, to ensure compa;bility. DON’T CHANGE THIS!! function [Yhistory, ErrorX, ErrorR, Iter]= Jacobi(L, a, YL, YR, M, tolX, tolR)

b) Create the A and b matrices that represent the M equa;ons for Yi in the form A Y = b where Y = [Y1 , … , YM ]T. Make sure you understand my example for M = 6 on the last page, before tackling this general M x M system..

c) For your ini;al guess, Y0 , just use the average value between the supports ½(YL+YR) for all Yi.

d) Use a “while” loop to keep itera;ng un;l the L2-norms of both the “proxy” error between itera;ons and the residual error are within their tolerances.

e) Calculate the column vector Yk for each itera;on k using the Jacobi itera;ve method.

f) Store each itera;on as a column vector in the variable Yhistory. So, for example, if you’re done a]er 1000 itera;ons, then Yhistory would be a 201 x 1000 matrix, with each column vector Yk represen;ng one itera;on. This will allow us to make plots of what the itera;ons looked like later.

g) Calculate (and store in vectors ErrorX and ErrorR) the L2-norms of both the “proxy” error in Yk and the residual error. This will allow us to make a convergence history plot later.
HW5 February 15

5.4 conTnued …

When it works, HW5_4 will automa;cally make two plots in one window. Save a .pdf of the plot window called Y.pdf. Then study the plots and your calculated variables and answer the following ques;ons:

i.  How many itera;ons did it take to converge? Hint: It took me many 10,000 s of itera;ons. ii.  For your final itera;on, what’s the lowest height of the cable (i.e minimum y(x))?

iii.  Your convergence history plot of “proxy error” should show an rapid ini;al improvement in error, and then a (rela;vely slow) linear convergence rate un;l comple;on. Use the slope in the linear region to approximate how many itera;ons it takes for the error to improve by each factor of 10 (i.e. gives you one more decimal of accuracy).

Submit online in the HW5 assignment:

•  Your working, documented Jacobi.m func;on, and plot Y.pdf.

•  Answers to ques;ons (i), (ii) and (iii) above in the online “comment” box.


HINTS for doing 5.4

1. 
Make sure you read and understand how I made my matrices for the M = 6 example. Then write

out on a piece of paper, before you even think of star;ng to use MATLAB, what your set of 201

equa;ons for 201 unknowns should generally look like. Keep the ordering of the equa;ons similar

to my example so your diagonal of A stays as “strong as possible”.
2. 
Sketch out on paper exactly what you need your Jacobi.m func;on to do! Break it up into high-

level jobs, then break that down into smaller parts, THEN start coding.
3. 
You should check each part of your code as you go along. For example, one thing you need to do is

create your A and b matrices. Don’t just jump into the M=201 system right away. Instead, try M=6

the first ;me and make sure you get exactly the same A and b matrices from my earlier example

using x = 2. If that’s correct, then it’s much more likely that your 201x201 system will work.
4. 
Is your actual Jacobi itera;on part working? (i.e. gegng Yk from the previous Yk-1) You should

double-check that part of your code alone against what you do know works. Take the 3x3 Ax=b

problem we did in class using Jacobi itera;on by-hand. We did the next 2 itera;ons on the board.

Your code should give exactly those results star;ng from those same A and b matrices.
5. 
What about your while loop? We already went over this a lot in HW3.7. I don’t fundamentally see

anything different in how you would set it up here, except we’re using norms of vectors to

represent the errors, rather than just scalar errors before.
6. 
Any confusion about norms? Well, try them out separately in MATLAB, or by hand. Give it a

smaller, known vector of errors as a test case and make sure you’re using the norm command for

L2-norms in MATLAB properly.
7. 
Does your code just keep running without stopping? Then get more informa;on WHY. Look at the

one of the later itera;ons – is it gegng closer to the final solu;on, or further away? Look at the

“proxy” and residual errors – are they gegng smaller, bigger, or kind of staying the same? Most problems here are caused by (a) you made incorrect A and b matrices and/or you did not carefully order the equa;ons to keep a strong diagonal (see hint 1 above), (b) your itera;on process is incorrect (see hint 4 above), (c) your while statement is not legng you end (see hint 5 above).
HW5 February 15

HINTS for doing 5.4 (conTnued …)

The boKom line is – these are the steps you need to take to develop your code properly the first ;me, and then to follow if you suspect something is not working. This is helping you be self-reliant, and not just reliant on the TAs. If you come to me or the TAs with a ques;on about the code not working, we’re happy to help, but come prepared! We’re going to ask you to show us that you already checked all the things above, and you’ll get turned away un;l you do so if you don’t bring us (for example) something on paper showing how you’re planning your code, or your output for trying to make the 6x6 A and b matrices from my earlier example, etc.

We do not want you to hand in all your preliminary papers or debugging work with your HW submission. That’s just for you to help you develop the code, and to communicate efficiently with us if there are any problems. Only hand in exactly what I asked for in problem 5.4.