Nume ical and Compu a ional Me hods II
(M´e odos Num´e icos y Compu acionales II)
Enginee ing Physics (UPC - ETSETB)
Compu a ional P ac icals
`
A. Mesegue
Enginye ia F´ısica Dep . F´ısica Aplicada (UPC)
Nume ical and Compu a ional Me hods II (2015-2016) A. Mesegue & F. Ma ques
Nume ical lab 1: Nume ical Linea Algeb a I (LU ac o iza ion).
1. The igu e below shows an inc easing se ies o ci cui s consis ing o iden ical esis o s (R= 3 Ω)
connec ed o a ba e y V= 1 V. The ci cui a he bo om is a gene aliza ion o he i s 3 ci cui s.
We wan o compu e he in ensi ies Iko each esis o as well as he equi alen esis ance Re
be ween poin s aand bo he n-ladde sys em o a bi a y n≥2.
11 2 123
123n
a
b
a
b
a
b
a
b
I1I3I5I2n−3
I2n−1
I2I4I6
I1I3I5I2n−3
I2n−2
(a) Ki cho : con ince you sel ha he op and bo om in ensi ies a each s ep ladde a e he
same (so we only ha e 2n−1 unkowns: I1, . . . , I2n−1), and ha :
2R I1+R I2=V
I1−I2−I3= 0
−R I2+ 2R I3+R I4= 0
I3−I4−I5= 0
.
.
..
.
..
.
..
.
.
−R I2n−4+ 2R I2n−3+R I2n−2= 0
I2n−3−I2n−2−I2n−1= 0
−R I2n−2+ 3R I2n−1= 0.
I we de ine I,V∈R2n−1as he ec o o unknown in ensi ies I≡[I1, I2, . . . , I2n−1]Tand
he igh hand side ec o associa ed wi h he ba e y V≡[V, 0,...,0]T, espec i ely, we
can ew i e he abo e sys em in ma ix o m:
2R R 0 0 0 0 ···
1−1−1 0 0 0 ···
0−R2R R 0 0 ···
0 0 1 −1−1 0 ···
.
.
..
.
. 0 .........
.
.
..
.
..
.
..
.
.−R2R R 0
.
.
..
.
..
.
..
.
..
.
. 1 −1−1
.
.
..
.
..
.
..
.
..
.
. 0 −R3R
I1
I2
I3
I4
I5
.
.
.
I2n−3
I2n−2
I2n−1
=
V
0
0
0
0
.
.
.
0
0
0
.
In Ma lab, he ma ix abo e can be gene a ed wi h he commands epma and spdiags:
B = [ epma ([-R;1],2*n-1,1), epma ([2*R;-1],2*n-1,1), epma ([R;-1],2*n-1,1)];
A = ull(spdiags(B,[-1:1],2*n-2,2*n-1)) ; A = [A ; [ze os(1,2*n-3) -R 3*R]];
(b) Using you own pa=lu code, sol e he sys em abo e o n= 10,20 and 30. Fo each case,
plo (k, Ik) and obse e ha he in ensi ies decay exponen ially ( y semilogy).
(c) Compu e Re(n) o he cases abo e and check ha l´ım
n→∞
Re(n) = R(1 + √3).
Enginye ia F´ısica Dep . F´ısica Aplicada (UPC)
Nume ical and Compu a ional Me hods II (2015-2016) A. Mesegue & F. Ma ques
Nume ical lab 2: Nume ical Linea Algeb a II (q ac o iza ion and he Leas Squa es P oblem).
The igu e on he igh shows he unc ion:
(x) = e−20(x+0.25)2+ 0.25 sin(30 x) e−20(x−0.25)2,
wi hin he domain [−1,1]. In mnc1, you s udied he
polynomial in e pola ion echnique o app oxima-
e unc ions by imposing an in e polan Πn (x)∈
Rn[x] o adop he alues o he unc ion a n+ 1
gi en nodes x0, x1,...,xn, i.e.:
Πn (xj) = (xj), j = 0,1,...,n.
x
y
-1 -0.5 0 0.5 1
-0.2
0
0.5
1
In his case you ha e o app oxima e (x) wi h polynomial expansions o a bi a y deg ee by aking
an also a bi a y (and g ea e ) numbe o sampling poin s.
(a) Conside he app oxima ion p(x)∈Rm−1[x] based on a linea combina ion o he i s m−1
monomials 1, x, x2,...,xm−1:
p(x) =
m−1
X
j=0
ajxj=a0+a1x1+a2x2+···+am−1xm−1,
and he alues j= (xj) o he unc ion (x) e alua ed a he se o nequispaced nodes in
[−1,1]: xj=−1+2j/(n−1), j = 0,1,...,n−1, wi h n > m. I we de ine a= [a0, a1,...,am−1]T
and = [ 0, 1,..., n−1]T, he p oblem has mcoe icien s o be de e mined and ncondi ions o
be sa is ied by he polynomial.
Since n > m, imposing p(xk) = (xk), k = 0,1,...,n−1, leads o an o e de e mined Vande monde-
ype n×msys em Va= whose explici o m is:
1x0x2
0··· xm−1
0
1x1x2
1··· xm−1
1
1x2x2
2··· xm−1
2
.
.
..
.
..
.
..
.
.
1xn−1x2
n−1··· xm−1
n−1
a0
a1
a2
.
.
.
am−1
=
0
1
2
.
.
.
n−1
The Vande monde ma ix abo e can be easily cons uc ed by using ande (x) Ma lab’s buil -
in command, whe e xis he ec o o nodes. Then apply he lipl ( lip le - igh ) on he
esul ing ma ix o eo de he columns and d op he m+ 1,...,n columns wi h V=V(:,1:m).
The ec o a ha minimizes he no m o he esidual k −Vakcan be compu ed by means
o he q ac o iza ion. Use you own q -code o sol e he Vande monde sys em o he cases
(n, m) = (14,7),(28,14),(28,20),(64,30). Fo each case, plo he leas squa es polynomial on a
ine g id. Obse e he p og essi e imp o emen in cap u ing he small oscilla ions on he igh .
You will also obse e a amilia phenomenon nea he ends o he in e al.
(b) Op ional: epea he compu a ion wi h Chebyche poin s xj= cos jπ
n−1, j = 0,1,...,n−1.
Enginye ia F´ısica Dep . F´ısica Aplicada (UPC)
Nume ical and Compu a ional Me hods II (2015-2016) A. Mesegue & F. Ma ques
Nume ical lab 3: Nume ical Linea Algeb a III (G am-Schmid eo hogonaliza ion).
In his expe imen , you will compa e h ee di e en me hods o o hogonaliza ion: classical G am-
Schmid (cgs), G am-Schmid wi h eo hogonaliza ion (gs ) and q . You will o hogonalize Hilbe
ma ices H(n) ha ha e he explici elemen s
Hij =1
i+j−1, i, j = 1,2, . . . , n.
A Hilbe ma ix is gene a ed wi h he Ma lab command hilb(n), whe e nis he dimension o
he ma ix. You need h ee ex e nal codes o o hogonalize he ma ix Hso hey p o ide he Q
o hogonalized columns o he ma ix. To measu e he ’quali y’ o he o hogonaliza ion, you will
compu e he 2-no m kI−QTQk2, whe e Iis he iden i y ma ix. I he o hogonaliza ion is well done,
his no m should be e y small. To do ha , you will make use o he no m command ( ype help no m
o see which o he no ms you can ask Ma lab o compu e).
(a) cgs: gene a e H(n) (wi h n= 9 o n= 15, o example) and moni o he inne p oduc s hqj,˜qki
o j≤k−1 o he esul ing no malized basis. A e e y s age k, check he no m o ˜qkand
obse e how i is dec easing. Do you obse e a de e io a ion o kI−QTQk2when you inc ease
n?.
(b) gs : same as in (a) bu now moni o he non-o hogonal oundo s sj=hqj,˜qki, j ≤k. Subs ac
hese p ojec ions as seen in he lec u es and ecompu e he q(2)
k ec o s. Re-check o hogonali y.
Do you obse e an imp o emen when compa ed wi h cgs?.
(c) q : epea wi h Ma lab’s q . Which o he h ee me hods pe o ms be e ?
(d) Op ional: compu e he condi ion numbe o he Hilbe ma ices using he cond(H) command
(again, ype help cond o check in which no m is his condi ion numbe compu ed). A geome-
ical in e p e a ion o he ill- condi ioning o a ma ix is based on he lack o o hogonali y o i s
columns ( ha gene a e he ank o he ep esen ed linea map). You can isualize his phenome-
non by: (1s ) no malizing each column o H, (2nd) compu ing HTH, i.e., (HTH)ij = cos(θij),
whe e θij is he ’angle’ be ween he iand jcolumns o H:
o j = 1:n
H(:,j) = H(:,j)/no m(H(:,j)); % No maliza ion
end
CS = lipud(H’*H) ;% CS con ains he cosine o he angles
con ou (CS,[0.9:.01:1.0]); colo ba ; % Con ou s o cos(angle) in [0.9,1.0]
Whe eas a well-condi ioned ma ix should ha e he con ou s concen a ed nea by he diagonal, an
ill-condi ioned one will show a sp ead o his con ou lines along all i s elemen s. Repea he plo wi h
he Qma ix compu ed wi h cgs,gs and q .
Enginye ia F´ısica Dep . F´ısica Aplicada (UPC)
Nume ical and Compu a ional Me hods II (2012-2013) A. Mesegue & F. Ma ques
Nume ical lab 4: Nume ical Linea Algeb a IV (K ylo i e a i e me hods–GMRES).
The igu e below shows he same ci cui as in he i s p ac ical: iden ical esis o s (R= 3 Ω) connec ed
o a ba e y V= 1 V. We wan o compu e he in ensi ies Ikup o a ixed ole ance and wi hou
explici ly cons uc ing he ma ix. As seen in he i s p ac ical, he in ensi ies we e dec easing om
le o igh in an exponen ial ashion. The e o e, we may no be in e es ed in he whole O(m)- ec o
o in ensi ies [I1I2··· I2m−1]T, bu jus in i s i s leading O(n)-componen s ( he “ ail ” o his ec o
will con ain Ij’s whose magni ude is below ou p esc ibed ole ance, hus negligible). Tha is, we a e
going o look o O(n)-low dimensional app oxima ions o an O(m) dimensional p oblem.
12 3 m
a
b
I1I3I5I2m−3
I2m−1
I2I4I6
I1I3I5I2m−3
I2m−2
(a) Ki cho : he solu ion ec o o in ensi ies I= [I1I2··· I2m−1]Thas o sa is y he 2m−1
equa ions
2R I1+R I2=V,
I1−I2−I3= 0,
2R I2j−3+R I2j−2−R I2j−4= 0,(j= 3,...,m),
I2j−3−I2j−2−I2j−1= 0,(j= 3,...,m),
3R I2m−1−RI2m−2= 0.
In he i s p ac ical, he le hand side o he equa ions abo e we e seen as a ma ix- ec o
p oduc . Fo he K ylo me hod, we only need he ac ion o he ma ix on he ec o : I−→ AI.
Gene a e a Ma lab unc ion ha o a gi en ec o I= [I1I2··· I2m−1]T, p o ides he le
hand side o he equa ions abo e. Example o m= 15:
unc ion Ax = A un(x)
R = 3 ; V = 1 ; m = 15; Ax = [] ;
Ax = [Ax ; 2*R*x(1) + R*x(2)]; Ax = [Ax ; x(1) - x(2) - x(3)] ;
o j = 3:m
Ax = [Ax ; 2*R*x(2*j-3) + R*x(2*j-2) - R*x(2*j-4)];
Ax = [Ax ; x(2*j-3) - x(2*j-2) - x(2*j-1)] ;
end
Ax = [Ax ; 3*R*x(2*m-1) - R*x(2*m-2)] ;
The ou pu o his unc ion is he ec o Ax. Check ha you ou pu is he same as he one
p o ided by he code o he i s p ac ical.
(b) Using you own gm es code, sol e he sys em abo e o ole ances kAx−bk ≤ 10−3, 10−4and
10−5. Check ha he i s componen s o he solu ion ec o coincide wi h he ones p o ided
by he ma ix e sion. Also check ha he numbe o A noldi i e a ions g ows as long as you
dec ease he ole ance.
No e: you gm es code has o ha e an inpu -ou pu syn ax like his:
[x,k] = mygm es(A un,b, ol,dimk yl),
whe e A un is he ac ion (in oked wi h @A un) unc ion abo e, bis he igh hand side o he sys em,
ol is he ole ance and dimk yl is he dimension o he K ylo space whe e you app oxima ion is
been compu ed, i.e., n, so in his case dimk yl ≤2m−1. Op ional: ype help gm es in Ma lab.
Enginye ia F´ısica Dep . F´ısica Aplicada (UPC)
Nume ical and Compu a ional Me hods II (2015-2016) A. Mesegue & F. Ma ques
Nume ical lab 5: Sys ems o Nonlinea Equa ions (New on me hod and con inua ion).
The igu e on he igh shows a double pendulum consis ing o wo small
sphe es o mass mconnec ed by massless igid ods o equal leng h `.
The sys em is o ced o o a e wi h a cons an angula speed ωa ound
a e ical axis passing ough he pi o p. The goal o his p ac ical is
o ge amilia wi h he concep o mul iplici y o solu ions a ising om
nonlinea physical sys ems. In pa icula , in his exe cise we wan o
ind he possible equilib ium solu ions o a gi en angula speed, i.e.,
con igu a ions whe e he angles φ1and φ2 emain cons an in ime.
Do i be o e s a ing he p ac ical !: apply New on’s laws and impose
equilib ium in a e e ence ame co- o a ing wi h he wo masses. By
di iding adial and e ical equa ions o each mass you should ge :
F1(φ1, φ2, α) = an(φ1)−α(2 sin φ1+ sin φ2)=0
F2(φ1, φ2, α) = an(φ2)−2α(sin φ1+ sin φ2)=0,
whe e α=`ω2/2g is a dimensionless pa ame e . Fo simplici y, we will
assume ha φ1∈[0, π/2) and φ2∈(−π/2, π/2).
m
m
ℓ
ℓ
φ1
φ2
~g
ω
p
(a) The con igu a ion (φ1, φ2) = (0,0) is always an equilib ium solu ion ∀α. Fo α∈[0,2], compu e
he de e minan o he jacobian de JF(φ1, φ2) o he sys em F= [F1F2]Ta (φ1, φ2) = (0,0).
Plo de JF(0,0) as a unc ion o α. Do you ind any indica ion ha new b anches o equilib ium
solu ions may eme ge o some alue(s) o α?. Compu e he jacobian ma ix ia ini e di e ences
(no analy ically) and use he Ma lab command de o compu e he de e minan .
(b) New on explo a ion: o ixed alues o α∈[0,2], s a you New on me hod om ini ial andom
seeds (φ(0)
1, φ(0)
2)∈[0, π/2) ×(−π/2, π/2) and check i he i e a ion con e ges o o he solu ions.
I ha is he case, s o e only hose con e ged solu ions which all wi hin he in e al (φ1, φ2)∈
[0, π/2) ×(−π/2, π/2). Plo he con e ged solu ions φ1(α) and φ2(α) in wo di e en igu es.
(c) Con inua ion (op ional): s a ing om any o he new solu ions ound, con inua e hem as a
unc ion o αand plo hem in wo di e en igu es. A wha alues o αa e hese solu ions
bo n?
Enginye ia F´ısica Dep . F´ısica Aplicada (UPC)
Nume ical and Compu a ional Me hods II (2013-2014) A. Mesegue & F. Ma ques
Nume ical lab 6: Gauss-Legend e app oxima ion.
In his p ac ical you a e asked o code an app oxima ion scheme o he Gauss-Legend e nodes. This
scheme has o p o ide: (1) gl nodes {x0,···, xn}, (2) quad a u e weigh s {w0,...,wn}, (3) ba ycen ic
weigh s {λ0,...,λn}and (4) di e en ia ion ma ix Dij o an a bi a y numbe N=n+ 1 o nodes.
Code a Ma lab unc ion [x,w,lamb,D] = myleg(N) ha p o ides he abo e elemen s o a gi en N.
(a) Nodes and weigh s: ollow he s eps o P oblem 3, using he companion ma ix echnique.
(b) Ba ycen ic weigh s: once you ha e compu ed {xj}and {wj} o j= 0,...,n, use he he magical
o mula o ob ain he ba ycen ic weigh s:
λk= (−1)kq(1 −x2
k)wk.
Tha is, i you sample an a bi a y unc ion (x) a he gl nodes, he in e pola ing polynomial
(in ba ycen ic o m) is:
Πn (x) =
n
X
j=0
λj j(x−xj)−1
n
X
j=0
λj(x−xj)−1
, o x6== xjand wi h j= (xj), j = 0,...,n.
(c) Di e en ia ion ma ix: ecall om
mnc1 ha , in gene al:
Dij =
λj
λi
1
xi−xj
(i6=j),
X
k6=j
1
xj−xk
(i=j).
Tes you ma ix by nume ically di e en ia ing Run-
ge’s unc ion (x) = (1 + x2)−1and compa ing wi h
i s exac de i a i e ′(x) on he in e al [−1,1]. I we
de ine D jas he nume ical de i a i e o (x) a node
xj, plo on a semilogy igu e xj e sus |D j− ′(xj)|.
Inc ease nand check ha he e is no Runge ins abi-
li y and ha he con e gence is geome ical.
Fo highe alues o n, es you di e en ia ion ma ix
on he unc ion you used in P ac ical 2:
(x) = e−20(x+0.25)2+ 0.25 sin(30 x) e−20(x−0.25)2,
whose analy ical de i a i e (in Ma lab o ma ) is:
d = -((40*(x + .25))./exp(20*(x+.25).^2))-...
(10*(x - 1/4).*sin(30*x))./exp(20*(x-1/4).^2)+...
((15/2)*cos(30*x))./exp(20*(x-1/4).^2);
(d) Op ional: de elop a simila ou ine o he Chebyche nodes. This is easie han be o e because
you ha e al eady s udied he ba ycen ic weigh s and di e en ia ion ma ix o hese nodes in
mnc1. Bu emembe ha he quad a u e weigh s a e only use ul o in eg als o he o m
Z1
−1
(x)
√1−x2dx.
You will need bo h ou ines o he incoming p ac icals.
Enginye ia F´ısica Dep . F´ısica Aplicada (UPC)
Nume ical and Compu a ional Me hods II (2016-2017) A. Mesegue & F. Ma ques
Nume ical lab 7: Ma ix Exponen ia ion, Eigen alues and Disc e e Fou ie T ans o m.
The igu e on he igh shows a sys em o mas-
ses connec ed ia linea sp ings. Le xi( ) and
˙xi≡ i( ) be he displacemen o he i- h
block om i s equilib ium posi ion and i s
k1k2k3k4
m1m2m3
x1x2x3
ins an aneous speed a ime , espec i ely. You can e i y ha :
˙x1
˙x2
˙x3
˙ 1
˙ 2
˙ 3
=
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
β11 β12 0 0 0 0
β21 β22 β23 000
0β32 β33 000
x1
x2
x3
1
2
3
,wi h
β11 =−(k1+k2)/m1
β12 =k2/m1, β21 =k2/m2
β22 =−(k2+k3)/m2
β23 =k3/m2, β32 =k3/m3
β33 =−(k3+k4)/m3.
Le z( ) = [x1x2x3 1 2 3]T∈R6be he s a e ec o a ime , and B he 6×6 ma ix o coe icien s
o he ode sys em abo e. Fo a gi en a bi a y ini ial condi ion z0=z(0), he solu ion is z( ) = eB z0.
Random ini ial posi ions xi(0) and eloci ies i(0) will in gene al lead o coupled oscilla ions wi h
some cha ac e is ic equencies. The pu pose o his p ac ical is o iden i y hese na u al equencies
wi h wo di e en me hods.
Nume ical alues: ake, o example:
mi= 1, k1= 1, k2=√2, k3=√3, k4= 4,z0= [0.281 0.033 −1.33 1.12 0.35 −0.299]T.
(a) Me hod I: acco ding o wha you s udied in mechanics, he no mal modes o equencies a e he
imagina y pa s o he spec um o eigen alues o B. Using Ma lab’s command eig, compu e
he eigen alues o B. Fo he alues abo e, you should ge :
ω1= 1.02236572722, ω2= 1.8826701636, ω3= 2.5889863102.
(b) Me hod II: o a sampling se o imes j=
j∆ , j = 0,1,2,...,N −1, you can compu e
zj=z( j) = eB jz0’exac ly’ by using Ma lab’s
command expm(B* j) a e e y j. Howe e , he
command expm is compu a ionally e y expensi-
e. You jus need o compu e eB∆ once:
zj= ej∆ Bz0= e∆ Be(j−1)∆ Bz0= e∆ Bzj−1.
In eg a e he sys em abo e wi h ∆ = 0.25 and
N= 400 and s o e x2( i) (o any o he o he
a iables). Plo x2( ) o isualize he oscilla ion.
Using you own d ou ine, compu e he e-
quency spec um o x2and ep esen i on a
semilogy plo .
You d should show clea peaks a he no mal equencies (eigen alues o me hod i) as seen
in he plo s. The e o e, by analyzing he oscilla ion o a physical sys em we can in e i s in e nal
p ope ies (he e he sp ing cons an s kiand masses mj). In p ac ice, his ype o ime equency
(o spa ial wa eleng h) analysis is used in many b anches o Physics in o de o iden i y he
inne s uc u e o ma e . O cou se, his p ac ical is a e y simpli ied model.
Enginye ia F´ısica Dep . F´ısica Aplicada (UPC)
Nume ical and Compu a ional Me hods II (2015-2016) A. Mesegue & F. Ma ques
Nume ical lab 8: Bounda y Value P oblem (Quan um Mechanical Squa e Well)
The igu e on he igh shows he squa e po en ial well:
V(x) = V0,|x|> a/2
0,|x|< a/2.
By in oducing he change o a iable u=x/a, he Sch ¨odin-
ge equa ion:
d2ψ
dx2=2m
~2[V(x)−E]ψ,
V(x) = V0
V= 0
a/2
−a/20x
E
can be w i en as:
−d2φ
du2+βΘ(|u| − 1/2)φ=βεφ,
wi h: φ(u)≡ψ(x(u)), β =2ma2V0
~2,ε≡E
V0
and Θ(z) = 1, z > 0
0, z < 0.
Nume ical alues: use he same se o pa ame e s as in Eisbe g & Resnick (see Appendix G), i.e.,
β= 64. Se L≈5 in you ans o ma ion map om (−1,1) −→ (−∞,+∞).
(a) Compu e he ene gies εno he h ee bounded eigens a es. Rep esen he associa ed φneigen-
unc ions.
(b) Compa e you esul s wi h he ones gi en by he analy ical solu ion ob ained in Eisbe g &
Resnick (Appendix H). To do ha , emembe ha he no malized ene gies Ende ined by e &
a e:
En= mEna2
2~2=1
2pεnβ,
and hey ha e o sa is y ei he one o he wo ascenden al equa ions ( ha you can sol e by
means o a simple New on me hod):
En an En= β
4− E2
no Enco En=− β
4− E2
n.
You will need a e y high numbe o poin s (N∼256 o mo e). This is due o he ac ha
he po en ial is un ealis ic and ma hema ically ill-de ined (discon inuous), he e o e we canno
expec geome ical con e gence.
